/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_PAM_s1_v2.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_PAM_s1_v2
{
debug deb;
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double t;
	double tPrim;
	double U;
	double V;
	double h;
	double epsilonF;
	double beta;

	unsigned int L;
};

struct variables {
	double lnF;
	//double ln;
	//double lmC;
	double lmF;
	double nF;			// ok
	double mu;			// ok
	double mF;			// ok
	double dF;			// ok

    double Energy;
	double magnetization;			// ok
 };

/*
 *********************************************************
 *********************** ANALOG **************************
 *********************************************************
*/

volatile double   n = 1.60;
volatile double   t = -1.;
volatile double   tPrim = 0.;
volatile double   U = 9.;
volatile double   V = -2.;
volatile double   h = 0.1;
volatile double   epsilonF = -4.;
volatile double   beta = 500;
volatile int      L = 256;

volatile double   lnF = -3.7948168;
volatile double   lmF = 0.3664242;
volatile double   ln = 0;
volatile double   lmC = 0;
volatile double   nF = 0.8749191;
volatile double   mu = -0.5637287;
volatile double   mF = 0.2149817;
volatile double   dF = 0.2032936;

volatile double def_n;
volatile double def_m;
volatile double Energy;

double fq(int s) {
   return (2.*Power(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
             Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))),2))/((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double fdqdnF(int s) {
   return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
             Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
           (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/
                 Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
             1.*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
             2.*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
         (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double fdqdmF(int s) {
   return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
             Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
           (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
             1.*s*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
             2.*s*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
         (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double fdqddF(int s) {
   return 0. + (4.*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
           Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
         (0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
           Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
           (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/
            (2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
       ((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double fE(double epsilonk, int s, int l) {
   return (epsilonk - 2*ln - lnF - 2*mu + epsilonF - (lmC + lmF + 2*h)*s -
           l*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)))/2.;
}

double fdEdlnF(double epsilonk, int s, int l) {
   return (-(l*(epsilonk + lnF - epsilonF - lmC*s + lmF*s)) -
           Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)))/
         (2.*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}

double fdEdlmC(double epsilonk, int s, int l) {
   return (s*(l*(epsilonk + lnF - epsilonF - lmC*s + lmF*s) -
             Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2))))/
         (2.*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}

double fdEdlmF(double epsilonk, int s, int l) {
   return -(s*(l*(epsilonk + lnF - epsilonF - lmC*s + lmF*s) +
           Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2))))/
      (2.*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}

double fdEdq(double epsilonk, int s, int l) {
   return -((l*Power(V,2))/Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}


/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.t = p.tPrim = p.U = p.V = p.h = p.epsilonF = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.lnF = v.lmF = v.nF = v.mu = v.mF = v.dF = 0;
	v.Energy = v.magnetization = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.lnF += q.lnF*d;
	v.lmF += q.lmF*d;
	v.nF += q.nF*d;
	v.mu += q.mu*d;
	v.mF += q.mF*d;
	v.dF += q.dF*d;
	v.Energy += q.Energy*d;
	v.magnetization += q.magnetization*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.lnF = q.lnF*d;
	v.lmF = q.lmF*d;
	v.nF = q.nF*d;
	v.mu = q.mu*d;
	v.mF = q.mF*d;
	v.dF = q.dF*d;
	v.Energy = q.Energy*d;
	v.magnetization = q.magnetization*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.t = q.t;
	v.tPrim = q.tPrim;
	v.U = q.U;
	v.V = q.V;
	v.h = q.h;
	v.epsilonF = q.epsilonF;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.t  p.tPrim  p.U  p.V  p.h  p.epsilonF  p.beta  p.L
	variables v; // variables  v.lnF  v.lnF  v.lmC  v.lmF  v.nF  v.mu  v.mC  v.mF  v.d  v.Energy v.magnetization
	double nF, dF, mF, mu, lnF, lmF; // contractions

	// Abbreviations - dependend on kx, ky
	double epsilonk;

public:
	double E[4];
	double dEdlnF[4];
	double dEdlmC[4];
	double dEdlmF[4];
	double dEdq[4];

	double q[4];
	double dqdnF[4];
	double dqdmF[4];
	double dqddF[4];

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fq(int s);
	double fdqdnF(int s);
	double fdqdmF(int s);
	double fdqddF(int s);

    // Energy
    double fE(double epsilonk, int s, int l);
    double fdEdlnF(double epsilonk, int s, int l);
    double fdEdlmF(double epsilonk, int s, int l);
    double fdEdq(double epsilonk, int s, int l);

	// REFRESH
	void refresh() {
		nF = v.nF;
		dF = v.dF;
		mF = v.mF;
		mu = v.mu;
		lnF = v.lnF;
		lmF = v.lmF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		for(int i=0; i<2; i++) {
			for(int j=0; j<2; j++) {
				q[j+2*i] 	 = fq(1-2*i);
				dqdnF[j+2*i] = fdqdnF(1-2*i);
				dqdmF[j+2*i] = fdqdmF(1-2*i);
				dqddF[j+2*i] = fdqddF(1-2*i);
			}
		}

		// zeby pozniej mozna bylo wypisac
		glob.qup =  q[0];
		glob.qdown = q[2];

	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		epsilonk = 2.*p.t*(cos(kx) + cos(ky)) + 4.*p.tPrim*cos(kx)*cos(ky);

		refresh();

		for(int i=0; i<2; i++) { // spins: up are E[0] E[1]
			for(int j=0; j<2; j++) { // 1,2 family
				E[j+2*i] 	  = fE(epsilonk, 1-2*i, 1-2*j);
				dEdlnF[j+2*i] = fdEdlnF(epsilonk, 1-2*i, 1-2*j);
				dEdlmF[j+2*i] = fdEdlmF(epsilonk, 1-2*i, 1-2*j);
				dEdq[j+2*i]   = fdEdq(epsilonk, 1-2*i, 1-2*j);
			}
		}
	}
};


/////////////////
///////////////// FullSimplify
/////////////////
/*
double Func::fq(int s) {
	return (1.414213562373095*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))))/Sqrt((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fdqdnF(int s) {
	return (1.414213562373095*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		     0.7071067811865475*(1 - 1.*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))))/Power((nF + mF*s)*(1 - 0.5*(nF + mF*s)),1.5);
}

double Func::fdqdmF(int s) {
	return (1.414213562373095*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		     0.7071067811865475*s*(1 - 1.*nF - 1.*mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))))/Power((nF + mF*s)*(1 - 0.5*(nF + mF*s)),1.5);
}

double Func::fdqddF(int s) {
	return 0. + (1.414213562373095*(0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) + Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/(2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/Sqrt((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fE(double epsilonk, int s, int l) {
	return (epsilonk - 2*ln - lnF - 2*mu + p.epsilonF - (lmC + lmF + 2*p.h)*s - l*Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/2.;
}

double Func::fdEdlnF(double epsilonk, int s, int l) {
	return (-(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s)) - Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/
			   (2.*Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

double Func::fdEdlmC(double epsilonk, int s, int l) {
	return (s*(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s) - Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2))))/
			   (2.*Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

double Func::fdEdlmF(double epsilonk, int s, int l) {
	return -(s*(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s) + Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2))))/
			   (2.*Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

double Func::fdEdq(double epsilonk, int s, int l) {
	return (-2*l*Power(p.V,2)*fq(s))/Sqrt(4*Power(p.V,2)*Power(fq(s),2) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2));
}
*/
double Func::fq(int s) {
	return (2.*Power(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))),2))/((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fdqdnF(int s) {
	return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
		     (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/
		           Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
		       1.*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		       2.*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
		   (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double Func::fdqdmF(int s) {
	return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
		     (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
		       1.*s*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		       2.*s*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
		   (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double Func::fdqddF(int s) {
	return 0. + (4.*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
	      (0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/
	         (2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
	    ((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fE(double epsilonk, int s, int l) {
	double ln=0, lmC=0;
	return (epsilonk - 2*ln - lnF - 2*mu + p.epsilonF - (lmC + lmF + 2*p.h)*s -
		     l*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/2.;
}

double Func::fdEdlnF(double epsilonk, int s, int l) {
	double lmC=0;
	return (-(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s)) -
		     Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/
		   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

double Func::fdEdlmF(double epsilonk, int s, int l) {
	double lmC=0;
	return -(s*(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s) +
	        Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2))))/
	   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}

double Func::fdEdq(double epsilonk, int s, int l) {
	double lmC=0;
	return -((l*Power(p.V,2))/Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
}


/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void setParVar(parameters &p, variables &v) {
		this->p = p;
		this->v = v;
		f->set(p, v);
	}

	//void run(variables &w);
	void run2(variables &w);
	void run_global(variables &w);

	double give_fq(int s) { return f->fq(s); }
	double give_fdqdnF(int s) { return f->fdqdnF(s); }
	double give_fdqdmF(int s) { return f->fdqdmF(s); }
	double give_fdqddF(int s) { return f->fdqddF(s); }
	double give_fE(double epsilonk, int s, int l) { return f->fE(epsilonk, s, l); }
	double give_fdEdlnF(double epsilonk, int s, int l) { return f->fdEdlnF(epsilonk, s, l); }
	double give_fdEdlmF(double epsilonk, int s, int l) { return f->fdEdlmF(epsilonk, s, l); }
	double give_fdEdq(double epsilonk, int s, int l) { return f->fdEdq(epsilonk, s, l); }

};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	//if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2.};

	for(int i=0; i<4; i++) {
		w.lnF += Etan[i]*f->dEdq[i]*f->dqdnF[i];
		w.lmF += Etan[i]*f->dEdq[i]*f->dqdmF[i];
		w.nF += Etan[i]*f->dEdlnF[i];
		w.mu += Etan[i];
		w.mF += Etan[i]*f->dEdlmF[i];
		w.dF += Etan[i]*f->dEdq[i]*f->dqddF[i];

		if(i == 0 || i == 1) { w.magnetization += Etan[i]; }
		else { w.magnetization -= Etan[i]; }

        //if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        //if(f->E[i]<0) { w.Energy += f->E[i]; }

		if(DEEP_DEBUG) { printf("########### %f, %f; %d: %f | %e, %e\n", x, y, i, Etan[i], f->E[i], p.beta*f->E[i]); }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E[i] > 30.) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
		else if(p.beta*f->E[i] < -30) { w.Energy += f->E[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E[i]))/p.beta; }
	}
	if(DEEP_DEBUG) {
		printf("###--- lnF=%.15e lmF%.15e nF=%.15e mu=%.15e mF=%.15e dF=%.15e || %.15e %.15e %.15e %.15e | %.15e\n",
				w.lnF, w.lmF, w.nF, w.mu, w.mF, w.dF, f->E[0], p.beta*f->E[0], -exp(-p.beta*f->E[0])/p.beta, -log(1 + exp(-p.beta*f->E[0]))/p.beta, Etan[0]);
	}
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //

void Equa::run2(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	int max = p.L/2; // attention: L should be even!

	//int iii = 0;

/*
    // SERIAL CASE - OK
    for(int i=-max; i<max; i++) {
		for(int j=-max; j<max; j++) {
			onePoint(temp, i, j);
     	}
    }
    addTo(temp, w, 1);
    clear(temp);
*/
    // SMART WAY - ALSO OK
    // SERIAL CASE
   	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
   	// * edges (4 times)
   	for(int i=1; i<max; i++) {
   		onePoint(temp, i, i);
   		onePoint(temp, max, i);
   		onePoint(temp, i, 0);
   		//iii += 12;
   	}
      addTo(temp, w, 4);
      clear(temp);
      // * interior area (8 times) - SERIAL CASE
      for(int i=1; i<max; i++) {
   		for(int j=1; j<i; j++) {
   			onePoint(temp, i, j); //iii+=8;									//	sleep(1);
        	}
   	}
   	addTo(temp, w, 8);
   	clear(temp);
   	// * missing corners
   	onePoint(w, max, max);
   	onePoint(w, max, 0);
   	onePoint(w, 0, max);
   	onePoint(w, 0, 0);
   	//iii += 4;

    //printf("%d\n", iii);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.lnF /= L2;
	w.lmF /= L2;
	w.nF /= L2;
	w.mu /= L2;
	w.mF /= L2;
	w.dF /= L2;
    w.Energy /= L2;
    w.magnetization /= L2;

	w.lnF += v.lnF;
	w.lmF += v.lmF;
	w.nF += v.nF;
	w.mu += -p.n;
	w.mF += v.mF;
	w.dF += 2.*p.U*v.dF;

    v.Energy = w.Energy + v.lnF*v.nF + v.lmF*v.mF + p.U*v.dF*v.dF;
    v.magnetization = w.magnetization;

    // DEBUG
    /*
    double kx = 0.3;
    double ky = 0.2;
    double s = 1;
    double l = 1;
    double epsilonk = 2.*p.t*(cos(kx) + cos(ky)) + 4.*p.tPrim*cos(kx)*cos(ky);
    deb.fq = give_fq(s);
    deb.fdqdnF = give_fdqdnF(s);
    deb.fdqdmF = give_fdqdmF(s);
    deb.fdqddF = give_fdqddF(s);
    deb.fE = give_fE(epsilonk, s, l);
    deb.fdEdlnF = give_fdEdlnF(epsilonk, s, l);
    deb.fdEdlmF = give_fdEdlmF(epsilonk, s, l);
    deb.fdEdq = give_fdEdq(epsilonk, s, l);

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }
*/

	if(DEBUG2) {
		printf("###+++ lnF=%.15e lmF%.15e nF=%.15e mu=%.15e mF=%.15e dF=%.15e\n",
				w.lnF, w.lmF, w.nF, w.mu, w.mF, w.dF);
	}

    // results
    case_PAM_s1_v2::glob.Energy = v.Energy + v.mu*p.n;
    case_PAM_s1_v2::glob.magnetization = v.magnetization;
	//glob.qt = f->fgt();
}

void Equa::run_global(variables &w)
{
	clear(w);

	// REFRESH

	case_PAM_s1_v2::n = p.n;
	case_PAM_s1_v2::t = p.t;
	case_PAM_s1_v2::tPrim = p.tPrim;
	case_PAM_s1_v2::U = p.U;
	case_PAM_s1_v2::V = p.V;
	case_PAM_s1_v2::h = p.h;
	case_PAM_s1_v2::epsilonF = p.epsilonF;
	case_PAM_s1_v2::beta = p.beta;
	case_PAM_s1_v2::L = p.L;

	case_PAM_s1_v2::lnF = v.lnF;
	case_PAM_s1_v2::lmF = v.lmF;
	case_PAM_s1_v2::nF = v.nF;
	case_PAM_s1_v2::mu = v.mu;
	case_PAM_s1_v2::mF = v.mF;
	case_PAM_s1_v2::dF = v.dF;

	case_PAM_s1_v2::def_n = 0;
	case_PAM_s1_v2::def_m = 0;
	case_PAM_s1_v2::Energy = 0;

	double y0_d=0, y1_nF=0, y2_lnF=0, y3_mu=0, y4_lmF=0, y5_mF=0;

	  for(int i=0; i<L; i++) {
	      for(int j=0; j<L; j++) {
	         for(int s=-1; s<=1; s+=2) {
	            for(int l=-1; l<=1; l+=2) {
	               double kx = (double)i*2.*M_PI/L;
	               double ky = (double)j*2.*M_PI/L;
	               double epsilonk = 2.*t*(cos(kx) + cos(ky)) + 4.*tPrim*cos(kx)*cos(ky);
	               double E = fE(epsilonk, s, l);
	               double fermi = (1-tanh(beta*E/2.))/2.;

	               y0_d +=     fermi*fdEdq(epsilonk, s, l)*fdqddF(s);
	               y1_nF +=    fermi*fdEdlnF(epsilonk, s, l);
	               y2_lnF +=   fermi*fdEdq(epsilonk, s, l)*fdqdnF(s);
	               y3_mu +=    fermi*(-1);
	               y4_lmF +=   fermi*fdEdq(epsilonk, s, l)*fdqdmF(s);
	               y5_mF +=    fermi*fdEdlmF(epsilonk, s, l);

	               def_n += fermi;
	               def_m += fermi*s;

	               if(beta*E > 30.) { Energy += -exp(-beta*E)/beta; } // basically 0
	               if(beta*E < -30) { Energy += E; }
	               else { Energy += -log(1 + exp(-beta*E))/beta; }

	               //printf("###%4d %4d %3d %3d: %e\n", i, j, s, l, Energy);
	            }
	         }
	      }
	  }
	  //printf("%f ---- ", Energy);

	  y0_d /=     L*L;
	  y1_nF /=    L*L;
	  y2_lnF /=   L*L;
	  y3_mu /=    L*L;
	  y4_lmF /=   L*L;
	  y5_mF /=    L*L;

	  def_n /= L*L;
	  def_m /= L*L;
	  Energy /= L*L;

	  y0_d +=     2.*U*dF;
	  y1_nF +=    nF;
	  y2_lnF +=   lnF;
	  y3_mu +=    n;
	  y4_lmF +=   lmF;
	  y5_mF +=    mF;

	  Energy += lnF*nF + lmF*mF + U*dF*dF + mu*n;

    // results
	w.lnF = y2_lnF;
	w.lmF = y4_lmF;
	w.nF = y1_nF;
	w.mu = y3_mu;
	w.mF = y5_mF;
	w.dF = y0_d;


    case_PAM_s1_v2::glob.Energy = Energy;
    case_PAM_s1_v2::glob.magnetization = def_m;
	//glob.qt = f->fgt();
    //printf("%f\n", Energy);

    return;
}

int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.t = ((Par *) params)->t();
	p.tPrim = ((Par *) params)->tPrim();
	p.U = ((Par *) params)->U();
	p.V = ((Par *) params)->V();
	p.h = ((Par *) params)->h();
	p.epsilonF = ((Par *) params)->epsilonF();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.lnF = gsl_vector_get (x, 0);
	v.lmF = gsl_vector_get (x, 1);
	v.nF = gsl_vector_get (x, 2);
	v.mu = gsl_vector_get (x, 3);
	v.mF = gsl_vector_get (x, 4);
	v.dF = gsl_vector_get (x, 5);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    //E.run(w);
    //E.run_global(w);
    E.run2(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.lnF);
	gsl_vector_set (f, 1, w.lmF);
	gsl_vector_set (f, 2, w.nF);
	gsl_vector_set (f, 3, w.mu);
	gsl_vector_set (f, 4, w.mF);
	gsl_vector_set (f, 5, w.dF);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# | E = %24.15e m = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			case_PAM_s1_v2::glob.Energy, case_PAM_s1_v2::glob.magnetization);
		printf("\n");
	}

    return GSL_SUCCESS;
}

void ismaximum(Par &inpar, Var &invar) {
	printf("#"); inpar.print_verboten();
	Var initV = invar;

	// Solving
	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var exact(*res.get_V());
	printf("#Wynik: \n#");
	res.print_label();
	root.print();
	printf("\n");
	initV = exact;

	// EXTRA
	//Var exact = invar;
	//initV = exact;

	// Checking
	printf("#is it maximum?\n");
	printf("#          lnF                  lmF                  nF                       mu                      mF                     d\n");

	struct parameters p;
	p.n = inpar.n();
	p.t = inpar.t();
	p.tPrim = inpar.tPrim();
	p.U = inpar.U();
	p.V = inpar.V();
	p.h = inpar.h();
	p.epsilonF = inpar.epsilonF();
	p.beta = inpar.beta();
	p.L = inpar.L();

	variables v; clear(v);
	Equa E(p, v);

	// EXTRA
	//v.lnF = initV.lnF();
	//v.lmF = initV.lmF();
	//v.nF = initV.nF();
	//v.mu = initV.mu();
	//v.mF = initV.mF();
	//v.dF = initV.dF();
	//variables w; //wynik
	//E.setParVar(p,v);
	//E.run2(w);

	double OptE = glob.Energy;
	double OptEnomu = glob.Energy - inpar.n()*initV.mu();
	while(1) {
		fflush(stdout);
		initV = exact; initV.shake(1e-2);
		initV.print();
		v.lnF = initV.lnF();
		v.lmF = initV.lmF();
		v.nF = initV.nF();
		v.mu = initV.mu();
		v.mF = initV.mF();
		v.dF = initV.dF();

		variables w; //wynik
		E.setParVar(p,v);
		E.run2(w);

		double globE = glob.Energy;
		double globEnomu = glob.Energy - p.n*v.mu;

		double dif = globE-OptE;
		double difnomu = globEnomu-OptEnomu;
		printf("%22.15f %22.15f | %13e", glob.Energy, glob.magnetization, dif); if(dif<0) { printf(" M!"); } else { printf("   "); }
		printf(" | (nomu): %13e", difnomu); if(difnomu<0) { printf(" M!"); }
		printf("\n");

		// EXTRA
		if(difnomu < 0) {
			OptE = globE; OptEnomu = globEnomu;
			initV.set_lnF(v.lnF);
			initV.set_lmF(v.lmF);
			initV.set_nF(v.nF);
			initV.set_mu(v.mu);
			initV.set_mF(v.mF);
			initV.set_dF(v.dF);
		}
	}

	printf("#end ismaximum\n\n\n");
	fflush(stdout);
	return;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, t=%f, tPrim=%f, U=%f, V=%f, h=%f, epsilonF=%f, beta=%f, L=%d\n", n(), t(), tPrim(), U(), V(), h(), epsilonF(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   t   tPrim   U   V   h   epsilonF  beta   L");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "            lnF                    lmF                    nF                       mu                      mF                     d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_PAM_s1_v2::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%22.15f %22.15f %22.15f %22.15f", glob.Energy, glob.magnetization, glob.qup, glob.qdown); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%22.15f %22.15f %22.15f\n", glob.Energy, glob.qup, glob.qdown); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                       E                       m                       qup                     qdown\n");}

void example() {

	printf("\nEXAMPLE PAM_s2\n");

    Par initPar_skosyB_HM(N_PAR);
    initPar_skosyB_HM.set_n(1.8);
    initPar_skosyB_HM.set_t(-1.);
    initPar_skosyB_HM.set_tPrim(0.25);
    initPar_skosyB_HM.set_U(5.);
    initPar_skosyB_HM.set_V(-0.5);
    initPar_skosyB_HM.set_h(0.0);
    initPar_skosyB_HM.set_epsilonF(-3.0);
    initPar_skosyB_HM.set_beta(5000.);
    initPar_skosyB_HM.set_L(512);
    //1.800000 -0.500000     -2.351769728783318e+00   2.011597641752943e+00   9.931289820265322e-01  -9.453773643233782e-01   9.724992321618903e-01   9.124046650986967e-02
    //-4.564498229273310      0.921813720205289      0.936344125127788      0.894677547172757
    Var initV_skosyB_HM(N_VAR);
    initV_skosyB_HM.set_lnF(-2.351769728783318e+00);
    initV_skosyB_HM.set_lmF(2.011597641752943e+00);
    initV_skosyB_HM.set_nF(9.931289820265322e-01);
    initV_skosyB_HM.set_mu(-9.453773643233782e-01);
    initV_skosyB_HM.set_mF(9.724992321618903e-01);
    initV_skosyB_HM.set_dF(9.124046650986967e-02);

    Par initPar_test2(N_PAR);
    initPar_test2.set_n(1.8);
    initPar_test2.set_t(-1.);
    initPar_test2.set_tPrim(0.25);
    initPar_test2.set_U(5.);
    initPar_test2.set_V(-1.2);
    initPar_test2.set_h(0.0);
    initPar_test2.set_epsilonF(-3.0);
    initPar_test2.set_beta(5000.);
    initPar_test2.set_L(512);
    //1.800000 -1.200000     -2.233440327819374e+00   2.425345289348903e-01   9.578033665563203e-01  -8.775002023231642e-01   2.574549076870326e-01   2.945306348908389e-01
    //-5.060748757804389      0.200000000000002      0.706205519159673      0.692179271627075
    Var initV_test2(N_VAR);
    initV_test2.set_lnF(-2.233440327819374e+00);
    initV_test2.set_lmF(2.425345289348903e-01);
    initV_test2.set_nF(9.578033665563203e-01);
    initV_test2.set_mu(-8.775002023231642e-01);
    initV_test2.set_mF(2.574549076870326e-01);
    initV_test2.set_dF(2.945306348908389e-01);



	// I'm choosing start point

	Par initPar(N_PAR); initPar = initPar_test2;  // PMh01;
	Var initV(N_VAR); initV = initV_test2; // PMh01;

	printf("#dochodze do punktu startowego\n");

	// STANDARD

	// FM up
	//automatic_run_epsilonF(initPar, initV, -2.64, -3.0, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 1.6, 1.80, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 8, 5, 0.5); printf("\n");
	//automatic_run_V(initPar, initV, -1.52, -1.2, 0.01); printf("\n");
	automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.01); printf("\n");


    printf("Koniec\n");
}

void example_inf() {
    Par initPar_mw(N_PAR);
    initPar_mw.set_n(1.);
    initPar_mw.set_t(-1.);
    initPar_mw.set_tPrim(0.25);
    initPar_mw.set_U(10.);
    initPar_mw.set_V(-2.);
    initPar_mw.set_h(0.0);
    initPar_mw.set_epsilonF(-3.);
    initPar_mw.set_beta(1500.);
    initPar_mw.set_L(512);
    initPar_mw.print_verboten();

	// I'm choosing start point
	Var initV(N_VAR);

	//FILE * pFile;
	//pFile = fopen("test2.txt", "w");

	// INFINITE
	Res res(N_VAR);
	RootGSL root(initPar_mw, function, res);
	printf("sdelta      "); res.print_label(); fflush(stdout);
	while(1) {
	    //initV.set_lnF(-3.118174067109039e+00 + 0.2*frand()-0.1);
	    //initV.set_lmF(8.498668693903373e-01 + 0.2*frand()-0.1);
	    //initV.set_nF(8.807308972838500e-01 + 0.2*frand()-0.1);
	    //initV.set_mu(-1.389603785358793e+00 + 1*frand()-0.5);
	    //initV.set_mF(4.651345420094439e-01 + 0.2*frand()-0.1);
	   // initV.set_dF(2.224796595854245e-01 + 0.2*frand()-0.1);

	    initV.set_lnF(18.*frand()-9.);
	    initV.set_lmF(10.*frand()-5.);
	    initV.set_nF(frand());
	    initV.set_mu(18.*frand()-9);
	    initV.set_mF(frand());
	    initV.set_dF(frand());


		if(root.Solve(initV) == GSL_SUCCESS) {
			//printf("       "); root.print(); root.print(pFile); par.print_verboten(); fflush(pFile); fflush(stdout);
			printf("       "); root.print();  fflush(stdout);

		}
	}
	//fclose (pFile);
}



void run_2D() {
	printf("#PAM s1 v2 - run 2D 7.01.2015\n");

    // (A) U=20; skosy; HM
    Par initPar_skosyA_HM(N_PAR);
    initPar_skosyA_HM.set_n(1.6);
    initPar_skosyA_HM.set_t(-1.);
    initPar_skosyA_HM.set_tPrim(0.25);
    initPar_skosyA_HM.set_U(20.);
    initPar_skosyA_HM.set_V(-1.76);
    initPar_skosyA_HM.set_h(0.0);
    initPar_skosyA_HM.set_epsilonF(-2.9);
    initPar_skosyA_HM.set_beta(1500.);
    initPar_skosyA_HM.set_L(512);
    //-2.900000 -1.760000     -2.420415528735513e+00   7.958824580537447e-01   7.655914641117010e-01  -1.053271023359923e+00
    //4.891182394422017e-01   7.847180487364411e-02    -5.103419306411800      0.399999999999976      0.737022209970833      0.483920411852615
    Var initV_skosyA_HM(N_VAR);
    initV_skosyA_HM.set_lnF(-2.420415528735513e+00);
    initV_skosyA_HM.set_lmF(7.958824580537447e-01);
    initV_skosyA_HM.set_nF(7.655914641117010e-01);
    initV_skosyA_HM.set_mu(-1.053271023359923e+00);
    initV_skosyA_HM.set_mF(4.891182394422017e-01);
    initV_skosyA_HM.set_dF(7.847180487364411e-02);

    // (A) U=20; skosy; SFM
    Par initPar_skosyA_FM(N_PAR);
    initPar_skosyA_FM.set_n(1.6);
    initPar_skosyA_FM.set_t(-1.);
    initPar_skosyA_FM.set_tPrim(0.25);
    initPar_skosyA_FM.set_U(20.);
    initPar_skosyA_FM.set_V(-0.5);
    initPar_skosyA_FM.set_h(0.0);
    initPar_skosyA_FM.set_epsilonF(-5.0);
    initPar_skosyA_FM.set_beta(1500.);
    initPar_skosyA_FM.set_L(512);
    //-5.000000 -0.500000     -6.793137287202764e+00   6.619212687443434e+00   9.937590235492925e-01
    //-1.292505478665554e+00   9.928243063840952e-01   1.653602049308368e-02    -6.303276526267961      0.980317474899311      0.976392500219461      0.663260295797182
    Var initV_skosyA_FM(N_VAR);
    initV_skosyA_FM.set_lnF(-6.793137287202764e+00);
    initV_skosyA_FM.set_lmF(6.619212687443434e+00);
    initV_skosyA_FM.set_nF(9.937590235492925e-01);
    initV_skosyA_FM.set_mu(-1.292505478665554e+00);
    initV_skosyA_FM.set_mF(9.928243063840952e-01);
    initV_skosyA_FM.set_dF(1.653602049308368e-02);

    // (A) U=20; skosy; PM
    Par initPar_skosyA_PM(N_PAR);
    initPar_skosyA_PM.set_n(1.6);
    initPar_skosyA_PM.set_t(-1.);
    initPar_skosyA_PM.set_tPrim(0.25);
    initPar_skosyA_PM.set_U(20.);
    initPar_skosyA_PM.set_V(-1.76);
    initPar_skosyA_PM.set_h(0.0);
    initPar_skosyA_PM.set_epsilonF(-2.90);
    initPar_skosyA_PM.set_beta(1500.);
    initPar_skosyA_PM.set_L(512);
    //1.600000 -1.760000 -2.900000    -2.295507050184221e+00  -4.564808146052940e-15   7.370023443414185e-01  -1.164486132841001e+00  -4.014837637427159e-15   8.172790716111851e-02    -5.049678734541228     -0.000000000000005      0.561667005163005      0.561667005163007    Var initV_skosyA_PM(N_VAR);
    Var initV_skosyA_PM(N_VAR);
    initV_skosyA_PM.set_lnF(-2.295507050184221e+00);
    initV_skosyA_PM.set_lmF(-4.564808146052940e-15);
    initV_skosyA_PM.set_nF(7.370023443414185e-01);
    initV_skosyA_PM.set_mu(-1.164486132841001e+00);
    initV_skosyA_PM.set_mF(-4.014837637427159e-15);
    initV_skosyA_PM.set_dF(8.172790716111851e-02);

    // (A) U=20; skosy; PM
    Par initPar_skosyA_PM35(N_PAR);
    initPar_skosyA_PM35.set_n(1.6);
    initPar_skosyA_PM35.set_t(-1.);
    initPar_skosyA_PM35.set_tPrim(0.25);
    initPar_skosyA_PM35.set_U(20.);
    initPar_skosyA_PM35.set_V(-3.5);
    initPar_skosyA_PM35.set_h(0.0);
    initPar_skosyA_PM35.set_epsilonF(-0.0);
    initPar_skosyA_PM35.set_beta(15000.);
    initPar_skosyA_PM35.set_L(512);
    //    0.000000 -3.500000     -2.099994297425928e+00  -4.459691661999458e-12   4.840152020044691e-01 -1.251194192117028e+00  -1.323882032593600e-12   1.076398953860862e-01    -5.590019372939784     -0.000000000001733      0.873656693459468      0.873656693459939
    Var initV_skosyA_PM35(N_VAR);
    initV_skosyA_PM35.set_lnF(-2.099994297425928e+00);
    initV_skosyA_PM35.set_lmF(-4.459691661999458e-12);
    initV_skosyA_PM35.set_nF(4.840152020044691e-01);
    initV_skosyA_PM35.set_mu(-1.251194192117028e+00);
    initV_skosyA_PM35.set_mF(-1.323882032593600e-12);
    initV_skosyA_PM35.set_dF(1.076398953860862e-01);

    // (B) U=8; skosy; HM
    Par initPar_skosyB_HM(N_PAR);
    initPar_skosyB_HM.set_n(1.6);
    initPar_skosyB_HM.set_t(-1.);
    initPar_skosyB_HM.set_tPrim(0.25);
    initPar_skosyB_HM.set_U(8.);
    initPar_skosyB_HM.set_V(-1.52);
    initPar_skosyB_HM.set_h(0.0);
    initPar_skosyB_HM.set_epsilonF(-2.64);
    initPar_skosyB_HM.set_beta(1500.);
    initPar_skosyB_HM.set_L(1024);
    //-2.640000 -1.520000     -1.972361479912702e+00   6.759122054833073e-01   8.039039560658447e-01  -1.224879890086142e+00
    //4.775435523450455e-01   1.652976399576606e-01    -4.786625786676127      0.399999999999992      0.807427684408615      0.675366174617013
    Var initV_skosyB_HM(N_VAR);
    initV_skosyB_HM.set_lnF(-1.972361479912702e+00);
    initV_skosyB_HM.set_lmF(6.759122054833073e-01);
    initV_skosyB_HM.set_nF(8.039039560658447e-01);
    initV_skosyB_HM.set_mu(-1.224879890086142e+00);
    initV_skosyB_HM.set_mF(4.775435523450455e-01);
    initV_skosyB_HM.set_dF(1.652976399576606e-01);

    // (B) U=8; skosy; SFM
    Par initPar_skosyB_FM(N_PAR);
    initPar_skosyB_FM.set_n(1.6);
    initPar_skosyB_FM.set_t(-1.);
    initPar_skosyB_FM.set_tPrim(0.25);
    initPar_skosyB_FM.set_U(8.);
    initPar_skosyB_FM.set_V(-0.5);
    initPar_skosyB_FM.set_h(0.0);
    initPar_skosyB_FM.set_epsilonF(-4.0);
    initPar_skosyB_FM.set_beta(1500.);
    initPar_skosyB_FM.set_L(1024);
    //-4.000000 -0.500000     -3.506264928239187e+00   3.275580122526352e+00   9.915246031089974e-01  -1.296788039426506e+00
    //9.857066869262093e-01   4.755205712310006e-02    -5.320265888126178      0.965210424861737      0.963205889652765      0.857228364106780
    Var initV_skosyB_FM(N_VAR);
    initV_skosyB_FM.set_lnF(-3.506264928239187e+00);
    initV_skosyB_FM.set_lmF(3.275580122526352e+00);
    initV_skosyB_FM.set_nF(9.915246031089974e-01);
    initV_skosyB_FM.set_mu(-1.296788039426506e+00);
    initV_skosyB_FM.set_mF(9.857066869262093e-01);
    initV_skosyB_FM.set_dF(4.755205712310006e-02);

    // (B) U=8; skosy; PM35
    Par initPar_skosyB_PM35(N_PAR);
    initPar_skosyB_PM35.set_n(1.6);
    initPar_skosyB_PM35.set_t(-1.);
    initPar_skosyB_PM35.set_tPrim(0.25);
    initPar_skosyB_PM35.set_U(8.);
    initPar_skosyB_PM35.set_V(-3.5);
    initPar_skosyB_PM35.set_h(0.0);
    initPar_skosyB_PM35.set_epsilonF(-0.0);
    initPar_skosyB_PM35.set_beta(1500.);
    initPar_skosyB_PM35.set_L(1024);
    //0.000000 -3.500000     -1.565502175652598e+00  -4.952900094883626e-14   5.438588619380501e-01
    //-1.662264817200630e+00  -1.539830531324181e-14   1.924397706082357e-01    -5.840259621285909     -0.000000000000011      0.949740114587498      0.949740114587500
    Var initV_skosyB_PM35(N_VAR);
    initV_skosyB_PM35.set_lnF(-1.565502175652598e+00);
    initV_skosyB_PM35.set_lmF(-4.952900094883626e-14);
    initV_skosyB_PM35.set_nF(5.438588619380501e-01);
    initV_skosyB_PM35.set_mu(-1.662264817200630e+00);
    initV_skosyB_PM35.set_mF(-1.539830531324181e-14);
    initV_skosyB_PM35.set_dF(1.924397706082357e-01);


    // (B) U=8; skosy; PM
    Par initPar_skosyB_PM(N_PAR);
    initPar_skosyB_PM.set_n(1.6);
    initPar_skosyB_PM.set_t(-1.);
    initPar_skosyB_PM.set_tPrim(0.25);
    initPar_skosyB_PM.set_U(8.);
    initPar_skosyB_PM.set_V(-1.52);
    initPar_skosyB_PM.set_h(0.0);
    initPar_skosyB_PM.set_epsilonF(-2.64);
    initPar_skosyB_PM.set_beta(15000.);
    initPar_skosyB_PM.set_L(1024);
    //1.600000 -1.520000 -2.640000    -1.857176298429466e+00   1.623103007003470e-12   7.780727250780457e-01  -1.276710501030139e+00   1.746121603420997e-12   1.730280278385541e-01    -4.747404470146168
    //1.600000 -1.520000 -2.640000    -1.857176298429466e+00   1.623103007003470e-12   7.780727250780457e-01  -1.276710501030139e+00   1.746121603420997e-12   1.730280278385541e-01    -4.747404470146168      0.000000000001761      0.688135070805229      0.688135070804720
    Var initV_skosyB_PM(N_VAR);
    initV_skosyB_PM.set_lnF(-1.857176298429466e+00);
    initV_skosyB_PM.set_lmF(1.623103007003470e-12);
    initV_skosyB_PM.set_nF(7.780727250780457e-01);
    initV_skosyB_PM.set_mu(-1.276710501030139e+00);
    initV_skosyB_PM.set_mF(1.746121603420997e-12);
    initV_skosyB_PM.set_dF(1.730280278385541e-01);


    // (C) U=20; skosy; HM
    Par initPar_skosyC_HM(N_PAR);
    initPar_skosyC_HM.set_n(1.6);
    initPar_skosyC_HM.set_t(-1.);
    initPar_skosyC_HM.set_tPrim(0.25);
    initPar_skosyC_HM.set_U(20.);
    initPar_skosyC_HM.set_V(-1.4);
    initPar_skosyC_HM.set_h(0.0);
    initPar_skosyC_HM.set_epsilonF(-1.6);
    initPar_skosyC_HM.set_beta(1500.);
    initPar_skosyC_HM.set_L(512);
//-1.600000 -1.400000     -1.439003263534609e+00   5.478879845754475e-01   6.787216972116955e-01  -8.171372384638584e-01
    //4.488986699991572e-01   5.486355566483660e-02    -3.772664176897585      0.399999999999994      0.804506614751801      0.527294921115259
    Var initV_skosyC_HM(N_VAR);
    initV_skosyC_HM.set_lnF(-1.439003263534609e+00);
    initV_skosyC_HM.set_lmF(5.478879845754475e-01);
    initV_skosyC_HM.set_nF(6.787216972116955e-01);
    initV_skosyC_HM.set_mu(-8.171372384638584e-01);
    initV_skosyC_HM.set_mF(4.488986699991572e-01);
    initV_skosyC_HM.set_dF(5.486355566483660e-02);

    // (C) U=20; skosy; SFM
    Par initPar_skosyC_FM(N_PAR);
    initPar_skosyC_FM.set_n(1.6);
    initPar_skosyC_FM.set_t(-1.);
    initPar_skosyC_FM.set_tPrim(0.25);
    initPar_skosyC_FM.set_U(20.);
    initPar_skosyC_FM.set_V(-0.5);
    initPar_skosyC_FM.set_h(0.0);
    initPar_skosyC_FM.set_epsilonF(-4.0);
    initPar_skosyC_FM.set_beta(1500.);
    initPar_skosyC_FM.set_L(512);
    //-4.000000 -0.500000     -5.331593536629276e+00   5.182142221836155e+00   9.897287253581268e-01
    //-1.279986222131751e+00   9.885857067735584e-01   1.574882317309905e-02    -5.311285975882894      0.974223323521100      0.975317479140524      0.536531077603267
    Var initV_skosyC_FM(N_VAR);
    initV_skosyC_FM.set_lnF(-5.331593536629276e+00);
    initV_skosyC_FM.set_lmF(5.182142221836155e+00);
    initV_skosyC_FM.set_nF(9.897287253581268e-01);
    initV_skosyC_FM.set_mu(-1.279986222131751e+00);
    initV_skosyC_FM.set_mF(9.885857067735584e-01);
    initV_skosyC_FM.set_dF(1.574882317309905e-02);

    // (C) U=20; skosy; PM
    Par initPar_skosyC_PM20(N_PAR);
    initPar_skosyC_PM20.set_n(1.6);
    initPar_skosyC_PM20.set_t(-1.);
    initPar_skosyC_PM20.set_tPrim(0.25);
    initPar_skosyC_PM20.set_U(20.);
    initPar_skosyC_PM20.set_V(-2.0);
    initPar_skosyC_PM20.set_h(0.0);
    initPar_skosyC_PM20.set_epsilonF(-0.0);
    initPar_skosyC_PM20.set_beta(1500.);
    initPar_skosyC_PM20.set_L(512);
    //0.000000 -2.000000     -1.199779568748108e+00  -1.143688361433124e-12   4.455996515322084e-01  -2.918217524339846e-01
    //-8.942477268874741e-13   6.501144605698526e-02    -3.628763682124841     -0.000000000000673      0.833135949338633      0.833135949339111
    Var initV_skosyC_PM20(N_VAR);
    initV_skosyC_PM20.set_lnF(-1.199779568748108e+00);
    initV_skosyC_PM20.set_lmF(-1.143688361433124e-12);
    initV_skosyC_PM20.set_nF(4.455996515322084e-01);
    initV_skosyC_PM20.set_mu(-2.918217524339846e-01);
    initV_skosyC_PM20.set_mF(-8.942477268874741e-13);
    initV_skosyC_PM20.set_dF(6.501144605698526e-02);

    Par initPar_skosyC_PM(N_PAR);
    initPar_skosyC_PM.set_n(1.6);
    initPar_skosyC_PM.set_t(-1.);
    initPar_skosyC_PM.set_tPrim(0.25);
    initPar_skosyC_PM.set_U(20.);
    initPar_skosyC_PM.set_V(-1.4);
    initPar_skosyC_PM.set_h(0.0);
    initPar_skosyC_PM.set_epsilonF(-1.6);
    initPar_skosyC_PM.set_beta(1500.);
    initPar_skosyC_PM.set_L(512);
    //1.600000 -1.400000 -1.600000    -1.357290810246480e+00   1.232858044121074e-13   6.367649490220507e-01  -7.132718675150935e-01   1.470621524107690e-13   5.724533374704660e-02    -3.753843752003188      0.000000000000139      0.637575258829808      0.637575258829719
    Var initV_skosyC_PM(N_VAR);
    initV_skosyC_PM.set_lnF(-1.357290810246480e+00);
    initV_skosyC_PM.set_lmF(1.232858044121074e-13);
    initV_skosyC_PM.set_nF(6.367649490220507e-01);
    initV_skosyC_PM.set_mu(-7.132718675150935e-01);
    initV_skosyC_PM.set_mF(1.470621524107690e-13);
    initV_skosyC_PM.set_dF(5.724533374704660e-02);

    // (D) U=20; skosy; HM
    Par initPar_skosyD_HM(N_PAR);
    initPar_skosyD_HM.set_n(1.6);
    initPar_skosyD_HM.set_t(-1.);
    initPar_skosyD_HM.set_tPrim(0.25);
    initPar_skosyD_HM.set_U(20.);
    initPar_skosyD_HM.set_V(-1.9);
    initPar_skosyD_HM.set_h(0.0);
    initPar_skosyD_HM.set_epsilonF(-3.6);
    initPar_skosyD_HM.set_beta(1500.);
    initPar_skosyD_HM.set_L(512);
    //-3.600000 -1.900000     -2.965300238348287e+00   9.040576339734441e-01   8.004420390795551e-01  -1.147528083816534e+00
    //5.058750405952265e-01   8.970624764277513e-02    -5.800911633616887      0.400000000000028      0.704277988163399      0.466573598877502
    Var initV_skosyD_HM(N_VAR);
    initV_skosyD_HM.set_lnF(-2.965300238348287e+00);
    initV_skosyD_HM.set_lmF(9.040576339734441e-01);
    initV_skosyD_HM.set_nF(8.004420390795551e-01);
    initV_skosyD_HM.set_mu(-1.147528083816534e+00);
    initV_skosyD_HM.set_mF(5.058750405952265e-01);
    initV_skosyD_HM.set_dF(8.970624764277513e-02);

    // (D) U=20; skosy; SFM
    Par initPar_skosyD_FM(N_PAR);
    initPar_skosyD_FM.set_n(1.6);
    initPar_skosyD_FM.set_t(-1.);
    initPar_skosyD_FM.set_tPrim(0.25);
    initPar_skosyD_FM.set_U(20.);
    initPar_skosyD_FM.set_V(-0.5);
    initPar_skosyD_FM.set_h(0.0);
    initPar_skosyD_FM.set_epsilonF(-6.0);
    initPar_skosyD_FM.set_beta(1500.);
    initPar_skosyD_FM.set_L(512);
    //-6.000000 -0.500000     -7.786840429545859e+00   7.580917231684396e+00   9.958375418403512e-01
    //-1.299913198481116e+00   9.949534092694577e-01   1.748637023519238e-02    -7.298173748303697      0.983309605077159      0.976051150090892      0.751589618042062
    Var initV_skosyD_FM(N_VAR);
    initV_skosyD_FM.set_lnF(-7.786840429545859e+00);
    initV_skosyD_FM.set_lmF(7.580917231684396e+00);
    initV_skosyD_FM.set_nF(9.958375418403512e-01);
    initV_skosyD_FM.set_mu(-1.299913198481116e+00);
    initV_skosyD_FM.set_mF(9.949534092694577e-01);
    initV_skosyD_FM.set_dF(1.748637023519238e-02);

    // (D) U=20; skosy; PM
    Par initPar_skosyD_PM40(N_PAR);
    initPar_skosyD_PM40.set_n(1.6);
    initPar_skosyD_PM40.set_t(-1.);
    initPar_skosyD_PM40.set_tPrim(0.25);
    initPar_skosyD_PM40.set_U(20.);
    initPar_skosyD_PM40.set_V(-4.0);
    initPar_skosyD_PM40.set_h(0.0);
    initPar_skosyD_PM40.set_epsilonF(-0.0);
    initPar_skosyD_PM40.set_beta(1500.);
    initPar_skosyD_PM40.set_L(512);
    //0.000000 -4.000000     -2.368730432073880e+00  -4.560837642420858e-13   4.951630856792805e-01
    //-1.612314883567463e+00  -1.209297134741038e-13   1.200064223946725e-01    -6.271885374472419     -0.000000000000168      0.884464548250999      0.884464548251037
    Var initV_skosyD_PM40(N_VAR);
    initV_skosyD_PM40.set_lnF(-2.368730432073880e+00);
    initV_skosyD_PM40.set_lmF(-4.560837642420858e-13);
    initV_skosyD_PM40.set_nF(4.951630856792805e-01);
    initV_skosyD_PM40.set_mu(-1.612314883567463e+00);
    initV_skosyD_PM40.set_mF(-1.209297134741038e-13);
    initV_skosyD_PM40.set_dF(1.200064223946725e-01);

    // (D) U=20; skosy; PM
    Par initPar_skosyD_PM(N_PAR);
    initPar_skosyD_PM.set_n(1.6);
    initPar_skosyD_PM.set_t(-1.);
    initPar_skosyD_PM.set_tPrim(0.25);
    initPar_skosyD_PM.set_U(20.);
    initPar_skosyD_PM.set_V(-1.9);
    initPar_skosyD_PM.set_h(0.0);
    initPar_skosyD_PM.set_epsilonF(-3.6);
    initPar_skosyD_PM.set_beta(1500.);
    initPar_skosyD_PM.set_L(512);
    //1.600000 -1.900000 -3.600000    -2.825369325012594e+00   7.679037764908704e-16   7.782565759817783e-01  -1.351079146645123e+00   7.983707112676729e-16   9.319992957224194e-02    -5.731764095470006      0.000000000000001      0.525901930629601      0.525901930629601
    Var initV_skosyD_PM(N_VAR);
    initV_skosyD_PM.set_lnF(-2.825369325012594e+00);
    initV_skosyD_PM.set_lmF(7.679037764908704e-16);
    initV_skosyD_PM.set_nF(7.782565759817783e-01);
    initV_skosyD_PM.set_mu(-1.351079146645123e+00);
    initV_skosyD_PM.set_mF(7.983707112676729e-16);
    initV_skosyD_PM.set_dF(9.319992957224194e-02);


    // (E) U=20; skosy; HM
    Par initPar_skosyE_HM(N_PAR);
    initPar_skosyE_HM.set_n(1.6);
    initPar_skosyE_HM.set_t(-1.);
    initPar_skosyE_HM.set_tPrim(0.25);
    initPar_skosyE_HM.set_U(20.);
    initPar_skosyE_HM.set_V(-1.2);
    initPar_skosyE_HM.set_h(0.0);
    initPar_skosyE_HM.set_epsilonF(-1.8);
    initPar_skosyE_HM.set_beta(15000.);
    initPar_skosyE_HM.set_L(512);
    //-1.800000 -1.200000     -1.351377622079572e+00   4.634197259496315e-01   7.166369306321255e-01  -7.859800242126275e-01   4.696414594263937e-01   4.825161237296344e-02    -3.704134954526555      0.400000000000130      0.757785611152762      0.460004203161496
    Var initV_skosyE_HM(N_VAR);
    initV_skosyE_HM.set_lnF(-1.351377622079572e+00);
    initV_skosyE_HM.set_lmF(4.634197259496315e-01);
    initV_skosyE_HM.set_nF(7.166369306321255e-01);
    initV_skosyE_HM.set_mu(-7.859800242126275e-01);
    initV_skosyE_HM.set_mF(4.696414594263937e-01);
    initV_skosyE_HM.set_dF(4.825161237296344e-02);

    // (E) U=20; skosy; SFM
    Par initPar_skosyE_FM(N_PAR);
    initPar_skosyE_FM.set_n(1.6);
    initPar_skosyE_FM.set_t(-1.);
    initPar_skosyE_FM.set_tPrim(0.25);
    initPar_skosyE_FM.set_U(20.);
    initPar_skosyE_FM.set_V(-0.5);
    initPar_skosyE_FM.set_h(0.0);
    initPar_skosyE_FM.set_epsilonF(-6.0);
    initPar_skosyE_FM.set_beta(15000.);
    initPar_skosyE_FM.set_L(512);
    //-6.000000 -0.500000     -7.786840429545859e+00   7.580917231684396e+00   9.958375418403512e-01
    //-1.299913198481116e+00   9.949534092694577e-01   1.748637023519238e-02    -7.298173748303697      0.983309605077159      0.976051150090892      0.751589618042062
    Var initV_skosyE_FM(N_VAR);
    initV_skosyE_FM.set_lnF(-7.786840429545859e+00);
    initV_skosyE_FM.set_lmF(7.580917231684396e+00);
    initV_skosyE_FM.set_nF(9.958375418403512e-01);
    initV_skosyE_FM.set_mu(-1.299913198481116e+00);
    initV_skosyE_FM.set_mF(9.949534092694577e-01);
    initV_skosyE_FM.set_dF(1.748637023519238e-02);

    // (E) U=20; skosy; PM
    Par initPar_skosyE_PM(N_PAR);
    initPar_skosyE_PM.set_n(1.6);
    initPar_skosyE_PM.set_t(-1.);
    initPar_skosyE_PM.set_tPrim(0.25);
    initPar_skosyE_PM.set_U(20.);
    initPar_skosyE_PM.set_V(-1.5);
    initPar_skosyE_PM.set_h(0.0);
    initPar_skosyE_PM.set_epsilonF(-0.0);
    initPar_skosyE_PM.set_beta(15000.);
    initPar_skosyE_PM.set_L(512);
    //0.000000 -1.500000     -8.735939855873387e-01   1.076807842158511e-12   4.317075236897177e-01  -4.821812339432058e-02   5.622655936026434e-13   4.869277037230212e-02    -3.018856811845896      0.000000000000674      0.815535735805682      0.815535735805334
    Var initV_skosyE_PM(N_VAR);
    initV_skosyE_PM.set_lnF(-8.735939855873387e-01);
    initV_skosyE_PM.set_lmF(1.076807842158511e-12);
    initV_skosyE_PM.set_nF(4.317075236897177e-01);
    initV_skosyE_PM.set_mu(-4.821812339432058e-02);
    initV_skosyE_PM.set_mF(5.622655936026434e-13);
    initV_skosyE_PM.set_dF(4.869277037230212e-02);


    // (F) U=20; skosy; HM
    Par initPar_skosyF_HM(N_PAR);
    initPar_skosyF_HM.set_n(1.6);
    initPar_skosyF_HM.set_t(-1.);
    initPar_skosyF_HM.set_tPrim(0.25);
    initPar_skosyF_HM.set_U(20.);
    initPar_skosyF_HM.set_V(-1.75);
    initPar_skosyF_HM.set_h(0.0);
    initPar_skosyF_HM.set_epsilonF(-3.0);
    initPar_skosyF_HM.set_beta(15000.);
    initPar_skosyF_HM.set_L(512);
    //-3.000000 -1.750000     -2.473211876462282e+00   7.970944301889716e-01   7.734603340326791e-01  -1.054582622367687e+00   4.936234573956655e-01   7.883688467292854e-02    -5.169708855452232      0.399999999999990      0.727387163173806      0.474065832391965
    Var initV_skosyF_HM(N_VAR);
    initV_skosyF_HM.set_lnF(-2.473211876462282e+00);
    initV_skosyF_HM.set_lmF(7.970944301889716e-01);
    initV_skosyF_HM.set_nF(7.734603340326791e-01);
    initV_skosyF_HM.set_mu(-1.054582622367687e+00);
    initV_skosyF_HM.set_mF(4.936234573956655e-01);
    initV_skosyF_HM.set_dF(7.883688467292854e-02);

    // (F) U=20; skosy; SFM
    Par initPar_skosyF_FM(N_PAR);
    initPar_skosyF_FM.set_n(1.6);
    initPar_skosyF_FM.set_t(-1.);
    initPar_skosyF_FM.set_tPrim(0.25);
    initPar_skosyF_FM.set_U(20.);
    initPar_skosyF_FM.set_V(-0.5);
    initPar_skosyF_FM.set_h(0.0);
    initPar_skosyF_FM.set_epsilonF(-6.0);
    initPar_skosyF_FM.set_beta(15000.);
    initPar_skosyF_FM.set_L(512);
    //-6.000000 -0.500000     -7.786840429545859e+00   7.580917231684396e+00   9.958375418403512e-01
    //-1.299913198481116e+00   9.949534092694577e-01   1.748637023519238e-02    -7.298173748303697      0.983309605077159      0.976051150090892      0.751589618042062
    Var initV_skosyF_FM(N_VAR);
    initV_skosyF_FM.set_lnF(-7.786840429545859e+00);
    initV_skosyF_FM.set_lmF(7.580917231684396e+00);
    initV_skosyF_FM.set_nF(9.958375418403512e-01);
    initV_skosyF_FM.set_mu(-1.299913198481116e+00);
    initV_skosyF_FM.set_mF(9.949534092694577e-01);
    initV_skosyF_FM.set_dF(1.748637023519238e-02);

    // (F) U=20; skosy; PM
    Par initPar_skosyF_PM30(N_PAR);
    initPar_skosyF_PM30.set_n(1.6);
    initPar_skosyF_PM30.set_t(-1.);
    initPar_skosyF_PM30.set_tPrim(0.25);
    initPar_skosyF_PM30.set_U(20.);
    initPar_skosyF_PM30.set_V(-3.0);
    initPar_skosyF_PM30.set_h(0.0);
    initPar_skosyF_PM30.set_epsilonF(-0.0);
    initPar_skosyF_PM30.set_beta(15000.);
    initPar_skosyF_PM30.set_L(512);
    //0.000000 -3.000000     -1.815086816502145e+00  -2.477154967833880e-12   4.720668290826348e-01  -9.067185953908885e-01  -8.862732797529284e-13   9.439882673257691e-02    -4.920199309844710     -0.000000000001046      0.861634577905505      0.861634577905864
    Var initV_skosyF_PM30(N_VAR);
    initV_skosyF_PM30.set_lnF(-1.815086816502145e+00);
    initV_skosyF_PM30.set_lmF(-2.477154967833880e-12);
    initV_skosyF_PM30.set_nF(4.720668290826348e-01);
    initV_skosyF_PM30.set_mu(-9.067185953908885e-01);
    initV_skosyF_PM30.set_mF(-8.862732797529284e-13);
    initV_skosyF_PM30.set_dF(9.439882673257691e-02);

    // (F) U=20; skosy; PM
    Par initPar_skosyF_PM(N_PAR);
    initPar_skosyF_PM.set_n(1.6);
    initPar_skosyF_PM.set_t(-1.);
    initPar_skosyF_PM.set_tPrim(0.25);
    initPar_skosyF_PM.set_U(20.);
    initPar_skosyF_PM.set_V(-1.75);
    initPar_skosyF_PM.set_h(0.0);
    initPar_skosyF_PM.set_epsilonF(-3.0);
    initPar_skosyF_PM.set_beta(15000.);
    initPar_skosyF_PM.set_L(512);
    //1.600000 -1.750000 -3.000000    -2.351110212834364e+00  -2.357829620383215e-13   7.467086244388889e-01  -1.186325047806285e+00  -3.883394075793493e-13   8.200473826980532e-02    -5.113006956663068     -0.000000000000298      0.549052609155344      0.549052609155534
    Var initV_skosyF_PM(N_VAR);
    initV_skosyF_PM.set_lnF(-2.351110212834364e+00);
    initV_skosyF_PM.set_lmF(-2.357829620383215e-13);
    initV_skosyF_PM.set_nF(7.467086244388889e-01);
    initV_skosyF_PM.set_mu(-1.186325047806285e+00);
    initV_skosyF_PM.set_mF(-8.862732797529284e-13);
    initV_skosyF_PM.set_dF(8.200473826980532e-02);

    // (F) U=20; skosy; PM
    Par initPar_skosyF_PMh01(N_PAR);
    initPar_skosyF_PMh01.set_n(1.6);
    initPar_skosyF_PMh01.set_t(-1.);
    initPar_skosyF_PMh01.set_tPrim(0.25);
    initPar_skosyF_PMh01.set_U(20.);
    initPar_skosyF_PMh01.set_V(-3.00);
    initPar_skosyF_PMh01.set_h(0.1);
    initPar_skosyF_PMh01.set_epsilonF(-0.0);
    initPar_skosyF_PMh01.set_beta(15000.);
    initPar_skosyF_PMh01.set_L(512);
    //0.100000 -3.000000 -0.000000    -1.819355855366787e+00   1.157030233955272e-01   4.732225532895803e-01  -9.084768687350150e-01   4.321841196328356e-02   9.434232361263011e-02    -4.923140112071909      0.058546496731693      0.870271860104055      0.852772627049093
    Var initV_skosyF_PMh01(N_VAR);
    initV_skosyF_PMh01.set_lnF(-1.819355855366787e+00);
    initV_skosyF_PMh01.set_lmF(1.157030233955272e-01);
    initV_skosyF_PMh01.set_nF(4.732225532895803e-01);
    initV_skosyF_PMh01.set_mu(-9.084768687350150e-01);
    initV_skosyF_PMh01.set_mF(4.321841196328356e-02);
    initV_skosyF_PMh01.set_dF(9.434232361263011e-02);

    // testy
    // (A) U=20; skosy; HM
    Par initPar_test(N_PAR);
    initPar_test.set_n(1.6);
    initPar_test.set_t(-1.);
    initPar_test.set_tPrim(0.25);
    initPar_test.set_U(5.);
    initPar_test.set_V(-1.6);
    initPar_test.set_h(0.0);
    initPar_test.set_epsilonF(-3.0);
    initPar_test.set_beta(15000.);
    initPar_test.set_L(512);
    //-3.000000 -1.600000
    Var initV_test(N_VAR);
    initV_test.set_lnF(-1.932415676210140e+00);
    initV_test.set_lmF(6.208234539206465e-01);
    initV_test.set_nF(8.873602617972758e-01);
    initV_test.set_mu(-1.679279295931211e+00);
    initV_test.set_mF(4.616483811185186e-01);
    initV_test.set_dF(2.732947996850061e-01);

	// I'm choosing start point

	Par initPar(N_PAR); initPar = initPar_skosyF_HM;  // PMh01;
	Var initV(N_VAR); initV = initV_skosyF_HM; // PMh01;

	printf("#dochodze do punktu startowego\n");

	// STANDARD

	// FM up
	//automatic_run_VtoP(initPar, initV, -0.5, -1.00, 0.01,  -12./5., -36./5.); printf("\n");
	automatic_run_n(initPar, initV, 1.6, 1.50, 0.01); printf("\n");
	//automatic_run_VtoP(initPar, initV, -1.0, -1.40, 0.01,  -12./5., -36./5.); printf("\n");

	// HM down
	///////automatic_run_VtoP(initPar, initV, -1.75, -1.20, 0.01,  -12./5., -36./5.); printf("\n");
	//automatic_run_n(initPar, initV, 1.6, 1.70, 0.01); printf("\n");
	//automatic_run_VtoP(initPar, initV, -1.75, -1.60, 0.01,  -12./5., -36./5.); printf("\n");

	// HM up & PM up
	//////automatic_run_VtoP(initPar, initV, -1.9, -1.2, 0.01,  -12./5., -36./5.); printf("\n");
	//automatic_run_n(initPar, initV, 1.60, 1.48, 0.005); printf("\n");
	//automatic_run_VtoP(initPar, initV, -1.75, -2.10, 0.01,  -12./5., -36./5.); printf("\n");

		//automatic_run_VtoP(initPar, initV, -0.5, -0.75, 0.01,  -12./5., -36./5.); printf("\n");
		//automatic_run_n(initPar, initV, 1.6, 1.57, 0.01); printf("\n");
		//automatic_run_U(initPar, initV, 20, 8.0, 1.0); printf("\n");
		//automatic_run_h(initPar, initV, 0., 0.1, 0.01); printf("\n");

		//automatic_run_V(initPar, initV, -0.5, -0.71, 0.01); printf("\n");
		//automatic_run_tPrim(initPar, initV, 0., 0.25, 0.05); printf("\n");
		//automatic_run_n(initPar, initV, 1.75, 2.50, 0.01); printf("\n");

		//automatic_run_epsilonF(initPar, initV, -3.0, -2.35, 0.05); printf("\n");


	//automatic_run_VtoP(initPar, initV, -3.0, -2.3, 0.01,  -12./5., -36./5.); printf("\n");

	printf("#Jestem w punkcie startowym\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // ADAPTACYJNY
    adaptic_2D_PMHM_nVtoP(initPar_skosyF_PM, initV_skosyF_PM, initV_skosyF_HM, 1.6, 2.0, 0.01, -1.75, -3.0, 0.001, -12./5., -36./5.);
    //adaptic_2D_PMHM_nVtoP(initPar_skosyF_PM, initV_skosyF_PM, initV_skosyF_HM, 1.6, 1.0, 0.001, -1.75, -3.0, 0.001, -12./5., -36./5.);

    //adaptic_2D_PMHM_nVtoP(initPar_skosyF_PM, initV_skosyF_PM, initV_skosyF_HM, 1.6, 2.0, 0.01, -1.75, -0.0, 0.01, -12./5., -36./5.);
    //adaptic_2D_PMHM_nVtoP(initPar_skosyF_PM, initV_skosyF_PM, initV_skosyF_HM, 1.6, 1.0, 0.01, -1.75, -0.0, 0.01, -12./5., -36./5.);


    // STANDARD
    //automatic_run_VtoP(initPar, initV, -1.40, -1.60, 0.0002,   -12./5., -36./5.); // FM up
    //automatic_run_VtoP(initPar, initV, -1.60, -1.40, 0.0002,   -12./5., -36./5.); // HM down
    //automatic_run_VtoP(initPar, initV, -2.10, -2.35, 0.0002,   -12./5., -36./5.); // HM up & PM up

    // Magnetyzacja (A)
    //automatic_run_VtoP(initPar, initV, -1.76, -0.00, 0.001, -5./3., -17.5/3); initPar=oldPar; initV=oldVar; automatic_run_VtoP(initPar, initV, -1.76, -3.50, 0.001, -5./3., -17.5/3);
    //automatic_run_VtoP(initPar, initV, -0.50, -3.50, 0.001, -5./3., -17.5/3);
    //automatic_run_VtoP(initPar, initV, -3.50, -0.00, 0.001, -5./3., -17.5/3);

    // Magnetyzacja (C)
    //automatic_run_VtoP(initPar, initV, -1.4, -0.00, 0.001, -8./3., -16./3.); initPar=oldPar; initV=oldVar; automatic_run_VtoP(initPar, initV, -1.4, -2.00, 0.001, -8./3., -16./3.);
    //automatic_run_VtoP(initPar, initV, -0.50, -2.00, 0.001, -8./3., -16./3.);
    //automatic_run_VtoP(initPar, initV, -2.00, -0.00, 0.001, -8./3., -16./3.);

    // Magnetyzacja (F)
    //automatic_run_VtoP(initPar, initV, -1.75, -0.00, 0.001, -12./5., -36./5.); initPar=oldPar; initV=oldVar; automatic_run_VtoP(initPar, initV, -1.76, -3.0, 0.001, -12./5., -36./5.);
    //automatic_run_VtoP(initPar, initV, -0.50, -3.0, 0.001, -12./5., -36./5.);
    //automatic_run_VtoP(initPar, initV, -3.0, -0.00, 0.001, -12./5., -36./5.);


    // SKOSY (A) HM
	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.76, -3.5, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.76, 0.0, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.76, -3.5, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;
	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.76, 0.0, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;

    // SKOSY (A) FM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -0.5, -2.0, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -0.5, -2.0, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;

    // SKOSY (A) PM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -3.5, -1.5, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -3.5, -1.5, 0.01, -5./3., -17.5/3); initPar=oldPar; initV=oldVar;


    // SKOSY (B) HM
	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.52, -3.5, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.52,  0.0, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.52, -3.5, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;
	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.52,  0.0, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;

    // SKOSY (B) FM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -0.5, -2.0, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -0.5, -2.0, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;

    // SKOSY (B) PM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -3.5, -1.3, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -3.5, -1.3, 0.01, -4./3., -14./3.); initPar=oldPar; initV=oldVar;


    // SKOSY (C) HM
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.4, -2.0, 0.01, -8./3., -16./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.4,  0.0, 0.01, -8./3., -16./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.4, -2.0, 0.01, -8./3., -16./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.4,  0.0, 0.01, -8./3., -16./3.); initPar=oldPar; initV=oldVar;

    // SKOSY (C) FM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -0.5, -2.0, 0.01, -8./3., -16./3.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -0.5, -2.0, 0.01, -8./3., -16./3.); initPar=oldPar; initV=oldVar;

    // SKOSY (C) PM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -3.5, -1.0, 0.01, -8./3., -16./3); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -3.5, -1.0, 0.01, -8./3., -16./3); initPar=oldPar; initV=oldVar;


    // SKOSY (D) HM
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.9, -4.0, 0.01, -12./7., -48./7.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.9,  0.0, 0.01, -12./7., -48./7.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.9, -4.0, 0.01, -12./7., -48./7.); initPar=oldPar; initV=oldVar;
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.9,  0.0, 0.01, -12./7., -48./7.); initPar=oldPar; initV=oldVar;

    // SKOSY (D) FM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -0.5, -2.0, 0.01, -12./7., -48./7.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -0.5, -2.0, 0.01, -12./7., -48./7.); initPar=oldPar; initV=oldVar;

    // SKOSY (D) PM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -4.0, -1.5, 0.01, -12./7., -48./7); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -4.0, -1.5, 0.01, -12./7., -48./7); initPar=oldPar; initV=oldVar;


    // SKOSY (F) HM
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 1.0, 0.01, -1.75, -2.3, 0.001, -12./5., -36./5.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 1.2, 0.01, -1.6,  -1.2, 0.0001, -12./5., -36./5.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.75, -3.0, 0.001, -12./5., -36./5.); initPar=oldPar; initV=oldVar;
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.02, -1.6,  -1.2, 0.0002, -12./5., -36./5.); initPar=oldPar; initV=oldVar;

    // SKOSY (F) FM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 1.2, 0.01, -1.2, -1.6, 0.0001, -12./5., -36./5.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.02, -1.2, -1.6, 0.0002, -12./5., -36./5.); initPar=oldPar; initV=oldVar;

    // SKOSY (F) PM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 1.0, 0.01, -2.3, -1.75, 0.001, -12./5., -36./5); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -3.0, -1.75, 0.001, -12./5., -36./5); initPar=oldPar; initV=oldVar;

	/////////////////////////////

    // SKOSY (E) HM
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.2, -1.5, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.2,  0.0, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.2, -1.5, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;
 	//automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.2,  0.0, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;

    // SKOSY (E) FM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -0.5, -1.5, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -0.5, -1.5, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;

    // SKOSY (E) PM
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 0.9, 0.01, -1.5, -0.5, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nVtoP(initPar, initV, 1.6, 2.0, 0.01, -1.5, -0.5, 0.01, -6., -9.); initPar=oldPar; initV=oldVar;

    // pole magnetyczne
    // 2D
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 0.0, 0.002, -1.75, -0.5, 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 0.0, 0.002, -1.75, -3., 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 1.0, 0.02, -1.75, -0.5, 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 1.0, 0.02, -1.75, -3., 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    // FM
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 0.0, 0.002, -0.75, -2., 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 1.0, 0.02, -0.75, -2., 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    // PM
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 0.0, 0.01, -3.00, -1.75, 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hVtoP(initPar, initV, 0.1, 1.0, 0.02, -3.00, -1.75, 0.002, -12./5., -36./5); initPar=oldPar; initV=oldVar; printf("\n");

    return;
}

// Przygotowanie do wykresu fazowego [h = 0] (1) [n=1.6, eF = -3] [uwaga! dla U=1000, mamy d na poziomie 0.001 ciągle!]
void run_1D_a() {


    // HALF-METAL
    Par initPar_0hm(N_PAR);
    initPar_0hm.set_n(1.6);
    initPar_0hm.set_t(-1.);
    initPar_0hm.set_tPrim(0.25);
    //initPar_hm.set_U(1e6);
    initPar_0hm.set_U(40);
    initPar_0hm.set_V(-2);
    initPar_0hm.set_h(0);
    initPar_0hm.set_epsilonF(-3.);
    initPar_0hm.set_beta(1500.);
    initPar_0hm.set_L(512);

    Var initV_0hm(N_VAR);
    initV_0hm.set_lnF(-2.663475727915975e+00);  // (-2.604159389773947e+00);
    initV_0hm.set_lmF(8.679034488962510e-01);  // (7.807520215643790e-01);
    initV_0hm.set_nF(7.429579316711511e-01);  // (7.346471391573127e-01);
    initV_0hm.set_mu(-1.030760087086375e+00);  // (-8.008005206766025e-01);
    initV_0hm.set_mF(4.838565431166796e-01);  // (4.888132921778925e-01);
    initV_0hm.set_dF(4.532034175810522e-02);  // (1.848574279628256e-03);

    // PARAMEGNET (h=0)
    Par initPar_0PM(N_PAR);
    initPar_0PM.set_n(1.6);
    initPar_0PM.set_t(-1.);
    initPar_0PM.set_tPrim(0.25);
    //initPar_0PM.set_U(1e6);
    initPar_0PM.set_U(40);
    initPar_0PM.set_V(-2);
    initPar_0PM.set_h(0.0);
    initPar_0PM.set_epsilonF(-3.);
    initPar_0PM.set_beta(1500.);
    initPar_0PM.set_L(512);

    Var initV_0PM(N_VAR);
    initV_0PM.set_lnF(-2.540533874884412e+00);	// (-2.553142919495844e+00);
    initV_0PM.set_lmF(-3.425808519291675e-15);	// (-1.937570963138407e-15);
    initV_0PM.set_nF(7.101493319909054e-01);	// (7.021425593706448e-01);
    initV_0PM.set_mu(-1.148430408315591e+00);	// (-1.058641598931175e+00);
    initV_0PM.set_mF(-1.559393638282477e-15);	// (-3.322552082410408e-15);
    initV_0PM.set_dF(4.729445240470302e-02);	// (1.808452020625865e-06);


	// I'm choosing start point
    Par initPar(N_PAR); initPar = initPar_0hm;
	Var initV(N_VAR); initV = initV_0hm;

	printf("#dochodze do punktu startowego\n");

	// STANDARD
	//automatic_run_h(initPar, initV, 0., 0.5, 0.01); printf("\n");
	//automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.05); printf("\n");
	automatic_run_V(initPar, initV, -2, -1.5, 0.1); printf("\n");
	automatic_run_U(initPar, initV, 40, 5.0, 1.); printf("\n");
	automatic_run_n(initPar, initV, 1.6, 1.3, 0.05); printf("\n");
	automatic_run_V(initPar, initV, -1.5, -1.05, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 1.2, 1.3, 0.05); printf("\n");

	//automatic_run_epsilonF(initPar, initV, -2, -4.0, 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 1., 0., 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 0.1, 0.01, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 1.75, 2.50, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 1, 12.0, 0.01); printf("\n");

	printf("#Jestem w punkcie startowym\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // STANDARD
    //automatic_run_n(initPar, initV, 1.6, 2.00, 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.1, 10., 0.01);
    //automatic_run_V(initPar, initV, -5, -5.1, 0.005); printf("\n"); initPar=oldPar; initV=oldVar;
    //automatic_run_V(initPar, initV, -2, -10.0, 0.1); printf("\n"); initPar=oldPar; initV=oldVar;

	//automatic_run_V(initPar, initV, -1.0, -0.5, 0.05); printf("\n");
	//automatic_run_V(initPar, initV, -0.5, -0.4, 0.001); printf("\n");
	//automatic_run_V(initPar, initV, -0.4, -0.0, 0.0001); printf("\n");

    // 2D
    automatic_2Drun_nV(initPar, initV, 1.3, 1.7, 0.01, -1.05, -0.85, 0.001); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nV(initPar, initV, 1.7, 1.9, 0.01, -2.3, -3.0, 0.001); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nV(initPar, initV, 1.6, 2.0, 0.01, -1.6, -5., 0.02); initPar=oldPar; initV=oldVar;
    //automatic_2Drun_nV(initPar, initV, 1.6, 0.5, 0.01, -1.6, -5., 0.02); initPar=oldPar; initV=oldVar;

    // 2D
    //automatic_2Drun_hV(initPar, initV, 0.1, 0.0, 0.01, -2, 0., 0.1); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.1, 2.0, 0.05, -2, 0., 0.1); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.1, 0.0, 0.01, -2, -10., 0.1); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.1, 2.0, 0.05, -2, -10., 0.1); initPar=oldPar; initV=oldVar; printf("\n");

	return;
}


// test dla  t=-1, t'=0.25, U=20, n=1.8, V=-1, epsilon=-2
// test dla marcina (od mojego punktu)
void run_1D_b() {


    // HALF-METAL
    Par initPar_0hm(N_PAR);
    initPar_0hm.set_n(1.6);
    initPar_0hm.set_t(-1.);
    initPar_0hm.set_tPrim(0.25);
    //initPar_hm.set_U(1e6);
    initPar_0hm.set_U(40);
    initPar_0hm.set_V(-2);
    initPar_0hm.set_h(0);
    initPar_0hm.set_epsilonF(-3.);
    initPar_0hm.set_beta(1500.);
    initPar_0hm.set_L(512);

    Var initV_0hm(N_VAR);
    initV_0hm.set_lnF(-2.663475727915975e+00);  // (-2.604159389773947e+00);
    initV_0hm.set_lmF(8.679034488962510e-01);  // (7.807520215643790e-01);
    initV_0hm.set_nF(7.429579316711511e-01);  // (7.346471391573127e-01);
    initV_0hm.set_mu(-1.030760087086375e+00);  // (-8.008005206766025e-01);
    initV_0hm.set_mF(4.838565431166796e-01);  // (4.888132921778925e-01);
    initV_0hm.set_dF(4.532034175810522e-02);  // (1.848574279628256e-03);

    // PARAMEGNET (h=0)
    Par initPar_0PM(N_PAR);
    initPar_0PM.set_n(1.6);
    initPar_0PM.set_t(-1.);
    initPar_0PM.set_tPrim(0.25);
    //initPar_0PM.set_U(1e6);
    initPar_0PM.set_U(40);
    initPar_0PM.set_V(-2);
    initPar_0PM.set_h(0.0);
    initPar_0PM.set_epsilonF(-3.);
    initPar_0PM.set_beta(1500.);
    initPar_0PM.set_L(512);

    Var initV_0PM(N_VAR);
    initV_0PM.set_lnF(-2.540533874884412e+00);	// (-2.553142919495844e+00);
    initV_0PM.set_lmF(-3.425808519291675e-15);	// (-1.937570963138407e-15);
    initV_0PM.set_nF(7.101493319909054e-01);	// (7.021425593706448e-01);
    initV_0PM.set_mu(-1.148430408315591e+00);	// (-1.058641598931175e+00);
    initV_0PM.set_mF(-1.559393638282477e-15);	// (-3.322552082410408e-15);
    initV_0PM.set_dF(4.729445240470302e-02);	// (1.808452020625865e-06);


	// I'm choosing start point
    Par initPar(N_PAR); initPar = initPar_0hm;
	Var initV(N_VAR); initV = initV_0hm;

	printf("#dochodze do punktu startowego\n");

	// STANDARD
	//automatic_run_h(initPar, initV, 0.1, 0.0, 0.01); printf("\n");
	//automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.05); printf("\n");
	automatic_run_n(initPar, initV, 1.6, 1.90, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 40, 80.0, 1); printf("\n");
	automatic_run_V(initPar, initV, -2, -3.0, 0.02); printf("\n");
	//automatic_run_epsilonF(initPar, initV, -2, -4.0, 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 1., 0., 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 0.1, 0.01, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 1.75, 2.50, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 1, 12.0, 0.01); printf("\n");
	automatic_run_L(initPar, initV, 512, 1024, 128); printf("\n");

	printf("#Jestem w punkcie startowym\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // STANDARD
    //automatic_run_n(initPar, initV, 1.6, 2.00, 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.1, 10., 0.01);
    automatic_run_V(initPar, initV, -3, -6.0, 0.01); printf("\n"); //initPar=oldPar; initV=oldVar;
    automatic_run_V(initPar, initV, -6, -5.0, 0.01); printf("\n"); initPar=oldPar; initV=oldVar;
    //automatic_run_V(initPar, initV, -2, -10.0, 0.1); printf("\n"); initPar=oldPar; initV=oldVar;

	return;
}

void run_1D_c() {


    // HALF-METAL
    Par initPar_0hm(N_PAR);
    initPar_0hm.set_n(1.6);
    initPar_0hm.set_t(-1.);
    initPar_0hm.set_tPrim(0.25);
    initPar_0hm.set_U(80);
    initPar_0hm.set_V(-6.4);
    initPar_0hm.set_h(0);
    initPar_0hm.set_epsilonF(-3.);
    initPar_0hm.set_beta(1500.);
    initPar_0hm.set_L(512);

    Var initV_0hm(N_VAR);
    initV_0hm.set_lnF(-5.989875438789134e+00);  // (-2.604159389773947e+00);
    initV_0hm.set_lmF(2.018320938540286e+00);  // (7.807520215643790e-01);
    initV_0hm.set_nF(6.003734094433115e-01);  // (7.346471391573127e-01);
    initV_0hm.set_mu(-3.160671739842151e+00);  // (-8.008005206766025e-01);
    initV_0hm.set_mF(3.181793805861438e-01);  // (4.888132921778925e-01);
    initV_0hm.set_dF(6.370881762880867e-02);  // (1.848574279628256e-03);

    Par initPar_0hm5(N_PAR);
    initPar_0hm5.set_n(1.6);
    initPar_0hm5.set_t(-1.);
    initPar_0hm5.set_tPrim(0.25);
    initPar_0hm5.set_U(80);
    initPar_0hm5.set_V(-5);
    initPar_0hm5.set_h(0);
    initPar_0hm5.set_epsilonF(-3.);
    initPar_0hm5.set_beta(1500.);
    initPar_0hm5.set_L(512);

    Var initV_0hm5(N_VAR);
    initV_0hm5.set_lnF(-4.967488788608874e+00);
    initV_0hm5.set_lmF(1.701648313149404e+00);
    initV_0hm5.set_nF(6.193650654850454e-01);
    initV_0hm5.set_mu(-2.371067655116553e+00);
    initV_0hm5.set_mF(3.496109858136962e-01);
    initV_0hm5.set_dF(5.135335120266389e-02);

    // PARAMEGNET (h=0)
    Par initPar_0PM(N_PAR);
    initPar_0PM.set_n(1.6);
    initPar_0PM.set_t(-1.);
    initPar_0PM.set_tPrim(0.25);
    //initPar_0PM.set_U(1e6);
    initPar_0PM.set_U(80);
    initPar_0PM.set_V(-6.4);
    initPar_0PM.set_h(0.0);
    initPar_0PM.set_epsilonF(-3.);
    initPar_0PM.set_beta(1500.);
    initPar_0PM.set_L(512);

    Var initV_0PM(N_VAR);
    initV_0PM.set_lnF(-5.583042897405084e+00);	//(-4.619565650003411e+00);
    initV_0PM.set_lmF(8.249007601035477e-14);	//(1.315960364859691e-14);
    initV_0PM.set_nF(5.647071308006496e-01);	//(5.792036149167114e-01);
    initV_0PM.set_mu(-3.211701842890607e+00);	//(-2.379344540205699e+00);
    initV_0PM.set_mF(1.404408387478755e-14);	//(2.818370861690415e-15);
    initV_0PM.set_dF(6.597609029024809e-02);	//(5.338234046254477e-02);

    Par initPar_0PM5(N_PAR);
    initPar_0PM5.set_n(1.6);
    initPar_0PM5.set_t(-1.);
    initPar_0PM5.set_tPrim(0.25);
    initPar_0PM5.set_U(80);
    initPar_0PM5.set_V(-5);
    initPar_0PM5.set_h(0.0);
    initPar_0PM5.set_epsilonF(-3.);
    initPar_0PM5.set_beta(1500.);
    initPar_0PM5.set_L(512);

    Var initV_0PM5(N_VAR);
    initV_0PM5.set_lnF(-4.619565650003413e+00);
    initV_0PM5.set_lmF(1.435456245847044e-15);
    initV_0PM5.set_nF(5.792036149167116e-01);
    initV_0PM5.set_mu(-2.379344540205696e+00);
    initV_0PM5.set_mF(3.802303375630077e-16);
    initV_0PM5.set_dF(5.338234046254482e-02);

    // WeakFM (h=0)
    Par initPar_0wFM(N_PAR);
    initPar_0wFM.set_n(1.6);
    initPar_0wFM.set_t(-1.);
    initPar_0wFM.set_tPrim(0.25);
    initPar_0wFM.set_U(80);
    initPar_0wFM.set_V(-5);
    initPar_0wFM.set_h(0.0);
    initPar_0wFM.set_epsilonF(-3.);
    initPar_0wFM.set_beta(1500.);
    initPar_0wFM.set_L(512);

    Var initV_0wFM(N_VAR);
    initV_0wFM.set_lnF(-4.792826906306350e+00);
    initV_0wFM.set_lmF(1.180127231815659e+00);
    initV_0wFM.set_nF(6.005324543539678e-01);
    initV_0wFM.set_mu(-2.392560204820462e+00);
    initV_0wFM.set_mF(2.557694779709815e-01);
    initV_0wFM.set_dF(5.232160054405983e-02);

    // U10 - granica PM
    Par initPar_0U10n13_FM(N_PAR);
    initPar_0U10n13_FM.set_n(1.3);
    initPar_0U10n13_FM.set_t(-1.);
    initPar_0U10n13_FM.set_tPrim(0.25);
    initPar_0U10n13_FM.set_U(10);
    initPar_0U10n13_FM.set_V(-2.3);
    initPar_0U10n13_FM.set_h(0.0);
    initPar_0U10n13_FM.set_epsilonF(-3.);
    initPar_0U10n13_FM.set_beta(1500.);
    initPar_0U10n13_FM.set_L(512);

    Var initV_0U10n13_FM(N_VAR);
    initV_0U10n13_FM.set_lnF(-2.199523441882246e+00);
    initV_0U10n13_FM.set_lmF(1.450383040751822e+00);
    initV_0U10n13_FM.set_nF(7.538385310759889e-01);
    initV_0U10n13_FM.set_mu(-2.574828461333933e+00);
    initV_0U10n13_FM.set_mF(6.217073482700813e-01);
    initV_0U10n13_FM.set_dF(1.251327102607874e-01);

    // U10 - granica PM
    Par initPar_0U10n13_PM(N_PAR);
    initPar_0U10n13_PM.set_n(1.3);
    initPar_0U10n13_PM.set_t(-1.);
    initPar_0U10n13_PM.set_tPrim(0.25);
    initPar_0U10n13_PM.set_U(10);
    initPar_0U10n13_PM.set_V(-2.3);
    initPar_0U10n13_PM.set_h(0.0);
    initPar_0U10n13_PM.set_epsilonF(-3.);
    initPar_0U10n13_PM.set_beta(1500.);
    initPar_0U10n13_PM.set_L(512);

    Var initV_0U10n13_PM(N_VAR);
    initV_0U10n13_PM.set_lnF(-1.702149081723531e+00);
    initV_0U10n13_PM.set_lmF(-1.244789418318668e-14);
    initV_0U10n13_PM.set_nF(6.505257901395891e-01);
    initV_0U10n13_PM.set_mu(-2.661039939459516e+00);
    initV_0U10n13_PM.set_mF(-7.287938603169225e-15);
    initV_0U10n13_PM.set_dF(1.494762987256232e-01);

    // U10 - granica PM
    Par initPar_0U10n07_FM(N_PAR);
    initPar_0U10n07_FM.set_n(0.7);
    initPar_0U10n07_FM.set_t(-1.);
    initPar_0U10n07_FM.set_tPrim(0.25);
    initPar_0U10n07_FM.set_U(10);
    initPar_0U10n07_FM.set_V(-6);
    initPar_0U10n07_FM.set_h(0.0);
    initPar_0U10n07_FM.set_epsilonF(-3.);
    initPar_0U10n07_FM.set_beta(1500.);
    initPar_0U10n07_FM.set_L(512);

    Var initV_0U10n07_FM(N_VAR);
    initV_0U10n07_FM.set_lnF(-1.053019600266540e+00);
    initV_0U10n07_FM.set_lmF(4.714596664077792e-01);
    initV_0U10n07_FM.set_nF(3.584268224163927e-01);
    initV_0U10n07_FM.set_mu(-7.501831957170138e+00);
    initV_0U10n07_FM.set_mF(1.754674517765249e-01);
    initV_0U10n07_FM.set_dF(1.006130943326988e-01);

    // U10 - granica PM
    Par initPar_0U10n07_PM(N_PAR);
    initPar_0U10n07_PM.set_n(0.7);
    initPar_0U10n07_PM.set_t(-1.);
    initPar_0U10n07_PM.set_tPrim(0.25);
    initPar_0U10n07_PM.set_U(10);
    initPar_0U10n07_PM.set_V(-6);
    initPar_0U10n07_PM.set_h(0.0);
    initPar_0U10n07_PM.set_epsilonF(-3.);
    initPar_0U10n07_PM.set_beta(1500.);
    initPar_0U10n07_PM.set_L(512);

    Var initV_0U10n07_PM(N_VAR);
    initV_0U10n07_PM.set_lnF(-1.033672608598541e+00);
    initV_0U10n07_PM.set_lmF(-1.171733504430959e-15);
    initV_0U10n07_PM.set_nF(3.495260661690555e-01);
    initV_0U10n07_PM.set_mu(-7.450672914705343e+00);
    initV_0U10n07_PM.set_mF(-6.345280414453025e-16);
    initV_0U10n07_PM.set_dF(1.115098712391216e-01);

    // WeakFM (1) (h=0)
    Par initPar_test_wFM(N_PAR);
    initPar_test_wFM.set_n(1.90);
    initPar_test_wFM.set_t(-1.);
    initPar_test_wFM.set_tPrim(0.25);
    initPar_test_wFM.set_U(40);
    initPar_test_wFM.set_V(-5.31);
    initPar_test_wFM.set_h(0.0);
    initPar_test_wFM.set_epsilonF(-3.);
    initPar_test_wFM.set_beta(1500.);
    initPar_test_wFM.set_L(1024);

    Var initV_test_wFM(N_VAR);
    initV_test_wFM.set_lnF(-5.963022878737513e+00); //(-5.613053663132222e+00); //(-5.542730287426636e+00); //(-5.542770217852462e+00); //(-5.542690729069892e+00); //(-4.213815694342416e+00); //(-6.743117484042414e+00); //(-6.900196391955999e+00); //(-5.277956319423748e+00); //(-5.132895845093229e+00);
    initV_test_wFM.set_lmF(5.503619587634548e-01); //(7.628098472338071e-01); //(-1.994181709122411e-02); //(-2.694481432462748e-02); //(-8.781525513670849e-03); //(-4.844179891453324e-03); //(2.284682787246277e-03); //(2.227988463170961e-03); //(-4.178993009364998e-03); //(-2.273007002375141e-03);
    initV_test_wFM.set_nF(6.564948567480838e-01); //(6.528843254887314e-01); //(6.465828790916689e-01); //(6.465865425291082e-01); //(6.465792503317731e-01); //(6.192855794593156e-01); //(6.625648858861791e-01); //(6.231395661797994e-01); //(6.382812016007640e-01); //(5.718360625504408e-01);
    initV_test_wFM.set_mu(-1.136616841882755e+00); //(-1.325302814434698e+00); //(-1.298860923920719e+00); //(-1.299248466648259e+00); //(-1.299037390032825e+00); //(-1.961788825284934e+00); //(-1.072930881530705e+00); //(-3.915606265882410e+00); //(-1.554633947049995e+00); //(-4.607146467944705e+00);
    initV_test_wFM.set_mF(1.057130363491234e-01); //(1.542476205457906e-01); //(-4.116102396144565e-03); //(-5.561493931098859e-03); //(-1.812579000966734e-03); //(-1.273553890883672e-03); //(3.829537411792709e-04); //(3.033851253645326e-04); //(-8.889175472174177e-04); //(-3.778592365251521e-04);
    initV_test_wFM.set_dF(1.272491542137446e-01); //(1.191438267346788e-01); //(1.197038193378008e-01); //(1.197034905587270e-01); //(1.197041362327919e-01); //(9.305552260670760e-02); //(1.453128568100947e-01); //(1.571386785792039e-01); //(1.149320532676589e-01); //(1.214923837969979e-01);

    //-5.310000   -5.963022878737513e+00   5.503619587634548e-01   6.564948567480838e-01
    //-1.136616841882755e+00   1.057130363491234e-01   1.272491542137446e-01   -10.068269809579714      0.099995302715350      0.773315390863489      0.732064774087751

    // WeakFM (1) (h=0)
    Par initPar_test_inne(N_PAR);
    initPar_test_inne.set_n(1.90);
    initPar_test_inne.set_t(-1.);
    initPar_test_inne.set_tPrim(0.25);
    initPar_test_inne.set_U(40);
    initPar_test_inne.set_V(-5.31);
    initPar_test_inne.set_h(0.0);
    initPar_test_inne.set_epsilonF(-3.);
    initPar_test_inne.set_beta(1500.);
    initPar_test_inne.set_L(1024);

    Var initV_test_inne(N_VAR);
    initV_test_inne.set_lnF(-5.928151151519982e+00); //(-5.277954252955670e+00); //(-5.132895318315496e+00);
    initV_test_inne.set_lmF(-2.394051691599596e-12); //(-5.494250356074972e-15); //(1.364301985605308e-11);
    initV_test_inne.set_nF(6.535870070636454e-01); //(6.382810011238570e-01); //(5.718360143971865e-01);
    initV_test_inne.set_mu(-1.124302638754816e+00); //(-1.554525469837018e+00); //(-4.607109596472323e+00);
    initV_test_inne.set_mF(-4.619614410394032e-13); //(-1.259729257974710e-15); //(2.267306332390689e-12);
    initV_test_inne.set_dF(1.275156542993931e-01); //(1.149320731361312e-01); //(1.214923906978095e-01);

    //-5.310000   -5.928151098418658e+00  -8.453769057582191e-13   6.535870026129653e-01
    //-1.124057676802678e+00  -1.612542524247148e-13   1.275156538990917e-01   -10.068233583420424     -0.000000000000146      0.751660627747495      0.751660627747557
    //1012   -5.928151151519982e+00  -2.394051691599596e-12   6.535870070636454e-01
    //-1.124302638754816e+00  -4.619614410394032e-13   1.275156542993931e-01   -10.068233504656394     -0.000000000000407      0.751660623639864      0.751660623640045

    // I'm choosing start point
    Par initPar(N_PAR); initPar = initPar_test_wFM;
	Var initV(N_VAR); initV = initV_test_wFM;
    //Par initPar(N_PAR); initPar = initPar_test_inne;
	//Var initV(N_VAR); initV = initV_test_inne;

	printf("#dochodze do punktu startowego\n");

	// STANDARD
	//automatic_run_h(initPar, initV, 0.1, 0.0, 0.01); printf("\n");
	//automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.05); printf("\n");
	//automatic_run_n(initPar, initV, 1.6, 1.90, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 40, 80.0, 1); printf("\n");
	//automatic_run_V(initPar, initV, -6.4, -5.0, 0.1); printf("\n");
	//automatic_run_epsilonF(initPar, initV, -2, -4.0, 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 1., 0., 0.1); printf("\n");
	//automatic_run_U(initPar, initV, 0.1, 0.01, 0.01); printf("\n");
	//automatic_run_n(initPar, initV, 1.75, 2.50, 0.01); printf("\n");
	//automatic_run_U(initPar, initV, 1, 12.0, 0.01); printf("\n");

	printf("#Jestem w punkcie startowym\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // STANDARD
    //automatic_run_n(initPar, initV, 1.6, 1.00, 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.1, 10., 0.01);
    automatic_run_V(initPar, initV, -5.31, -5.70, 0.01); printf("\n"); initPar=oldPar; initV=oldVar;
    //automatic_run_V(initPar, initV, -2, -10.0, 0.1); printf("\n"); initPar=oldPar; initV=oldVar;

    //automatic_run_L(initPar, initV, 512, 4800, 20);
    //automatic_run_L(initPar, initV, 512, 1024, 100);
    //automatic_run_L(initPar, initV, 1024, 8*1024, 1024);
    //automatic_run_L(initPar, initV, 8*1024, 32*1024, 4096);
	return;
}

// Przygotowanie do wykresu fazowego [h = 0] (1) [n=1.6, eF = -3] [uwaga! dla U=1000, mamy d na poziomie 0.001 ciągle!]
void run_1D_d() {

    // HALF-METAL
    Par initPar_03hm(N_PAR);
    initPar_03hm.set_n(1.6);
    initPar_03hm.set_t(-1.);
    initPar_03hm.set_tPrim(0.25);
    initPar_03hm.set_U(40);
    initPar_03hm.set_V(-3);
    initPar_03hm.set_h(0.03);
    initPar_03hm.set_epsilonF(-3.);
    initPar_03hm.set_beta(1500.);
    initPar_03hm.set_L(512);

    Var initV_03hm(N_VAR);
    initV_03hm.set_lnF(-3.443152231334297e+00);
    initV_03hm.set_lmF(1.219505802064852e+00);
    initV_03hm.set_nF(6.871496882722989e-01);
    initV_03hm.set_mu(-1.467329774919196e+00);
    initV_03hm.set_mF(4.233911171656275e-01);
    initV_03hm.set_dF(6.388269193168816e-02);
    //0.030000 -3.000000     -3.443152231334297e+00   1.219505802064852e+00   6.871496882722989e-01  -1.467329774919196e+00   4.233911171656275e-01   6.388269193168816e-02

    // PARAMEGNET
    Par initPar_03PM(N_PAR);
    initPar_03PM.set_n(1.6);
    initPar_03PM.set_t(-1.);
    initPar_03PM.set_tPrim(0.25);
    initPar_03PM.set_U(40);
    initPar_03PM.set_V(-6);
    initPar_03PM.set_h(0.03);
    initPar_03PM.set_epsilonF(-3.);
    initPar_03PM.set_beta(1500.);
    initPar_03PM.set_L(512);

    Var initV_03PM(N_VAR);
    initV_03PM.set_lnF(-5.063775564890757e+00);	// (-2.553142919495844e+00);
    initV_03PM.set_lmF(3.280950565731572e-01);	// (-1.937570963138407e-15);
    initV_03PM.set_nF(5.958264095819987e-01);	// (7.021425593706448e-01);
    initV_03PM.set_mu(-3.341649370702580e+00);	// (-1.058641598931175e+00);
    initV_03PM.set_mF(6.107921561201747e-02);	// (-3.322552082410408e-15);
    initV_03PM.set_dF(1.162443388718409e-01);	// (1.808452020625865e-06);
    //0.030000 -6.000000     -5.063775564890757e+00   3.280950565731572e-01   5.958264095819987e-01  -3.341649370702580e+00   6.107921561201747e-02   1.162443388718409e-01   -10.398709731548127      0.076562854142848      0.804554427721354      0.779923195717257

    // I'm choosing start point
    Par initPar(N_PAR); initPar = initPar_03hm;
	Var initV(N_VAR); initV = initV_03hm;

	printf("#dochodze do punktu startowego\n");

	// STANDARD
	//automatic_run_n(initPar, initV, 1.6, 1.40, 0.01);
	//automatic_run_V(initPar, initV, -3, -1.8, 0.1);
	//automatic_run_h(initPar, initV, 0.03, 0.01, 0.01);
	//automatic_run_h(initPar, initV, 0.03, 0.2, 0.01);
	//automatic_run_V(initPar, initV, -3, -2, 0.1);
	//automatic_run_U(initPar, initV, 40, 5, 1);
	//automatic_run_h(initPar, initV, 0.2, 0.1, 0.01);
	//automatic_run_V(initPar, initV, -3., -2.7, 0.05);
	//automatic_run_h(initPar, initV, 0.15, 0.1, 0.01);
	//automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.05); printf("\n");
	//automatic_run_V(initPar, initV, -3., -6., 0.01); printf("\n");
	//automatic_run_epsilonF(initPar, initV, -2, -4.0, 0.1); printf("\n");

	printf("#Jestem w punkcie startowym\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // 1D
/*    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.03, 0.05, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.03, 0.07, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.03, 0.09, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.03, 0.11, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.03, 0.13, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.03, 0.15, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.01, 0.02, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;

    automatic_run_h(initPar, initV, 0.01, 0.01, 0.01);
    automatic_run_V(initPar, initV, -2, -3., 0.001); initPar=oldPar; initV=oldVar;
*/
    //////

    //automatic_run_V(initPar, initV, -3, 0., 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_V(initPar, initV, -6, 0., 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.03, 0.2, 0.0005); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.002, 0.25, 0.002); initPar=oldPar; initV=oldVar;

    // 03HM
////    automatic_2Drun_hV(initPar, initV, 0.1, 0.0, 0.005, -2.7,  -2.2, 0.02); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.1, 0.3, 0.005, -1.8,  0., 0.02); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.005, -1.8, -5., 0.02); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.3, 0.005, -1.8, -5., 0.02); initPar=oldPar; initV=oldVar; printf("\n");

    // 03PM
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.005, -5, -1., 0.02); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.3, 0.005, -5, -1., 0.02); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.005, -5, -8., 0.1); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.5, 0.005, -5, -8., 0.1); initPar=oldPar; initV=oldVar; printf("\n");

    //automatic_2Drun_hV(initPar, initV, 0.03, 0.2, 0.001, -6, -3., 0.02); initPar=oldPar; initV=oldVar; printf("\n");

    automatic_run_V(initPar, initV, -3, -3.1, 0.01); initPar=oldPar; initV=oldVar;
    automatic_run_V(initPar, initV, -3, -2.9, 0.01); initPar=oldPar; initV=oldVar;
	return;
}


/////////////////////////////////////////////////////////////////////////
///*******************************************************************///
/////////////////////////////////////////////////////////////////////////

void automatic_run_n(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_n, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_n(start + ((double)i)*step);
		printf("%f ", inpar.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_n(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n\n");
	fflush(stdout);
	return;
}

void automatic_run_tPrim(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrim, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim       "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim(start + ((double)i)*step);
		printf("%f ", inpar.tPrim());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrim\n\n\n");
	fflush(stdout);

	inpar.set_tPrim(stop);
	invar = initV;
	return;
}

void automatic_run_U(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_U, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_U(start + ((double)i)*step);
		printf("%f ", inpar.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_U\n\n\n");
	fflush(stdout);

	inpar.set_U(stop);
	invar = initV;
	return;
}

void automatic_run_V(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_V, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_V(start + ((double)i)*step);
		printf("%f ", inpar.V());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_V\n\n\n");
	fflush(stdout);

	inpar.set_V(stop);
	invar = initV;
	return;
}

void automatic_run_h(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_h, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#h           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_h(start + ((double)i)*step);
		printf("%f ", inpar.h());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_h\n\n\n");
	fflush(stdout);

	inpar.set_h(stop);
	invar = initV;
	return;
}

void automatic_run_epsilonF(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_epsilonF, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#epsilonF    "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_epsilonF(start + ((double)i)*step);
		printf("%f ", inpar.epsilonF());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_epsilonF\n\n\n");
	fflush(stdout);

	inpar.set_epsilonF(stop);
	invar = initV;
	return;
}

void automatic_run_beta(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_beta, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#beta        "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_beta(start + ((double)i)*step);
		printf("%f ", inpar.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_beta\n\n\n");
	fflush(stdout);

	inpar.set_beta(stop);
	invar = initV;
	return;
}

void automatic_run_L(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_L, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#L           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_L(start + ((double)i)*step);
		printf("%d ", inpar.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_L\n\n\n");
	fflush(stdout);

	inpar.set_L(stop);
	invar = initV;
	return;
}


///2D
void automatic_2Drun_nV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_n(start1);
	par.set_V(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n      V        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_n(start1 + ((double)i)*step1);
		par.set_V(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			printf("%f %f   ", par.n(), par.V());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_n(stop1); par.set_V(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}


///2D
void automatic_2Drun_efV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_epsilonF(start1);
	par.set_V(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#eF     V        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_epsilonF(start1 + ((double)i)*step1);
		par.set_V(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			printf("%f %f   ", par.epsilonF(), par.V());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_epsilonF(stop1); par.set_V(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

///2D
void automatic_2Drun_Vef(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_V(start1);
	par.set_epsilonF(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#   V       eF        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_V(start1 + ((double)i)*step1);
		par.set_epsilonF(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_epsilonF(start2 + ((double)j)*step2);
			printf("%f %f   ", par.V(), par.epsilonF());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_V(stop1); par.set_epsilonF(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_Vn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_V(start1);
	par.set_n(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V      n        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_V(start1 + ((double)i)*step1);
		par.set_n(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_n(start2 + ((double)j)*step2);
			printf("%f %f   ", par.V(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_V(stop1); par.set_n(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_hV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_h(start1);
	par.set_V(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#h        V        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_h(start1 + ((double)i)*step1);
		par.set_V(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			printf("%f %f   ", par.h(), par.V());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_h(stop1); par.set_V(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_hVtoP(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double A, double B) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_h(start1);
	par.set_V(start2);
	par.set_epsilonF(A*start2+B);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#h          V          eF     "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_h(start1 + ((double)i)*step1);
		par.set_V(start2);
		par.set_epsilonF(A*start2+B);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			par.set_epsilonF(A*(start2 + ((double)j)*step2)+B);
			printf("%f %f %f  ", par.h(), par.V(), par.epsilonF());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_h(stop2); par.set_V(stop2); par.set_epsilonF(A*stop2+B);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

//cisnienie
void automatic_2Drun_VtoPn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double A, double B) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_V(start1);
	par.set_epsilonF(A*start1+B);
	par.set_n(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V          eF          n     "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_V(start1 + ((double)i)*step1);
		par.set_epsilonF(A*(start1 + ((double)i)*step1)+B);
		par.set_n(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_n(start2 + ((double)j)*step2);
			printf("%f %f %f  ", par.V(), par.epsilonF(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_V(stop1); par.set_epsilonF(A*stop1+B); par.set_n(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}


//cisnienie
void automatic_2Drun_nVtoP(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2, double A, double B) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_n(start1);
	par.set_V(start2);
	par.set_epsilonF(A*start2+B);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n          V          eF     "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_n(start1 + ((double)i)*step1);
		par.set_V(start2);
		par.set_epsilonF(A*start2+B);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			par.set_epsilonF(A*(start2 + ((double)j)*step2)+B);
			printf("%f %f %f  ", par.n(), par.V(), par.epsilonF());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); printf("\n");fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_n(stop2); par.set_V(stop2); par.set_epsilonF(A*stop2+B);
	printf("\n#END\n");
	fflush(stdout);

	return;
}



void automatic_run_VtoP(Par &inpar, Var &invar, double start, double stop, double step, double A, double B) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	inpar.set_V(start);
	inpar.set_epsilonF(A*start+B);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_V, start %f, stop %f, step %f, p(V) = %f*V + %f | ", start, stop, step, A, B); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V            eF     "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_V(start + ((double)i)*step);
		inpar.set_epsilonF(A*(start + ((double)i)*step)+B);
		printf("%f  %f", inpar.V(), inpar.epsilonF());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();printf("\n");
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_V\n\n\n");
	fflush(stdout);

	inpar.set_V(stop); inpar.set_epsilonF(A*stop+B);
	invar = initV;
	return;
}

// adaptacyjny
void adaptic_2D_PMHM_nVtoP(Par &par, Var &initV_PM, Var &initV_HM, double start_n, double stop_n, double step_n, double start_V, double stop_V, double step_V, double A, double B) {
	Res res_PM(N_VAR);
	RootGSL *root_PM = new RootGSL(par, function, res_PM);
	Res res_HM(N_VAR);
	RootGSL *root_HM = new RootGSL(par, function, res_HM);

	// set the beginning
	par.set_n(start_n);
	par.set_V(start_V);
	par.set_epsilonF(A*start_V+B);

	// print the head
	printf("#"); par.print_verboten();
	printf("#start_n=%f, stop_n=%f, step_n=%f, start_V=%f, stop_V=%f, step_V=%f\n", start_n, stop_n, step_n, start_V, stop_V, step_V);
	printf("#init_p(1)   "); initV_PM.print(); printf("\n");
	printf("#init_p(2)   "); initV_HM.print(); printf("\n");
	printf("#n          V          eF       step      "); res_PM.print_label(); //printf("\n");
	fflush(stdout);

	initV_PM.shake();
	initV_HM.shake();

	// I check the starting point and copying to the OLD variables
	int sucess_PM = root_PM->Solve(initV_PM);
	Var old_PM_lv1(*res_PM.get_V());
	int sucess_HM = root_HM->Solve(initV_HM);
	Var old_HM_lv1(*res_HM.get_V());
	if(sucess_PM != GSL_SUCCESS || sucess_HM != GSL_SUCCESS) { printf("# shot is not good.\n");  exit(EXIT_FAILURE); }

	int LOOPS = int((stop_n-start_n)/step_n);
	if(LOOPS < 0) { LOOPS = -LOOPS; step_n = -step_n; }
	for(int i=0; i<=LOOPS; i++) {
		double current_n = start_n + ((double)i)*step_n;
		par.set_n(current_n);
		par.set_V(start_V);
		par.set_epsilonF(A*start_V+B);

		// restoring the the variables from the old-container
		initV_PM = old_PM_lv1; initV_PM.shake();
		initV_HM = old_HM_lv1; initV_HM.shake();

		// solving and coping to old variables of velel 1 and 2
		int sucess_PM = root_PM->Solve(initV_PM);
			old_PM_lv1 = initV_PM;
			Var old_PM_lv2(*res_PM.get_V());
			double OLD_PM_E_lv2 = res_PM.glob.Energy;
			double OLD_PM_m_lv2 = res_PM.glob.magnetization;
		int sucess_HM = root_HM->Solve(initV_HM);
			old_HM_lv1 = initV_HM;
			Var old_HM_lv2(*res_HM.get_V());
			double OLD_HM_E_lv2 = res_HM.glob.Energy;
			double OLD_HM_m_lv2 = res_HM.glob.magnetization;
		if(sucess_PM != GSL_SUCCESS || sucess_HM != GSL_SUCCESS) {
			printf("# first shot for n = %f is not good (PM=%d, HM=%d, SUCESS=%d).\n", current_n, sucess_PM, sucess_PM, GSL_SUCCESS);
			initV_PM = old_PM_lv1;
			initV_HM = old_HM_lv1;
			continue;
		}

		// printing the first point in line
		printf("%f  %f  %f  %f", par.n(), par.V(), par.epsilonF(), step_V);
		root_PM->print();
		printf(" -1  ");
		root_HM->print(); printf("\n"); fflush(stdout);

		// the second loop
		if((stop_V < start_V && step_V > 0) || (stop_V > start_V && step_V < 0)) { step_V = -step_V; }
		double current_V = start_V;
		int stop_loop_lv2 = 1;
		if(start_V > stop_V) {stop_loop_lv2 = -1; }
		while(1) {
				if(current_V*stop_loop_lv2 >= stop_V*stop_loop_lv2) { printf("#end_V\n\n"); break; }

				par.set_V(current_V + step_V);
				par.set_epsilonF(A*(current_V+step_V)+B);

				initV_PM = old_PM_lv2; initV_PM.shake();
				initV_HM = old_HM_lv2; initV_HM.shake();

				int sucess_PM = root_PM->Solve(initV_PM);
				double PM_E = res_PM.glob.Energy;
				double PM_m = res_PM.glob.magnetization;
				int sucess_HM = root_HM->Solve(initV_HM);
				double HM_E = res_HM.glob.Energy;
				double HM_m = res_HM.glob.magnetization;

				// if there is problem to obtain solution, take smaller step
				if(sucess_PM != GSL_SUCCESS || sucess_HM != GSL_SUCCESS) {
					if(myabs(step_V/adapt_step_reduction) >=  adapt_min_step) {
						step_V /= adapt_step_reduction;
						printf("#NO GSL_SUCCESS -> smaller step and again! (n=%f, PM=%d, HM=%d, SUCESS=%d).\n", current_n, sucess_PM, sucess_PM, GSL_SUCCESS);
						continue;
					} else {
						printf("#NO GSL_SUCCESS -> step too small (%f, %e). Break the misery!\n\n", step_V, step_V/adapt_step_reduction);
						break;
					}
				}

				// if too big jump of magnetization, we might jump to other phase
				if(myabs(HM_m-OLD_HM_m_lv2) > adapt_max_m_difference) {
					if(myabs(step_V/adapt_step_reduction) >=  adapt_min_step) {
						step_V /= adapt_step_reduction;
						printf("#Big jump of magnetization -> smaller step and again! ");
						root_PM->print(); printf(" <PM|HM>  "); root_HM->print(); printf("\n"); fflush(stdout);
						continue;
					} else {
						printf("#Big jump of magnetization and the step too small (%f, %e). Break the misery!\n\n", step_V, step_V/adapt_step_reduction);
						break;
					}
				}

				// special case, if PM starts winning
				if(PM_E < HM_E && myabs(PM_m - HM_m) > adapt_min_m_difference) {
					if(myabs(step_V/adapt_step_reduction) >=  adapt_min_step) {
						step_V /= adapt_step_reduction;
						printf("#PM_E < HM_E -> smaller step! n=%f V=%f eF=%f step=%f: ", par.n(), par.V(), par.epsilonF(), step_V);
						root_PM->print(); printf(" <PM|HM>  "); root_HM->print(); printf("\n"); fflush(stdout);
						continue;
					} else {
						printf("#"); root_PM->print(); printf(" <PM|HM>  "); root_HM->print(); printf("\n"); fflush(stdout);
						printf("#PM_E < HM_E -> I got I order transition between HM and PM\n");
						printf("#n=%f, V in (%f, %f); eF in (%f, %f)\n\n", current_n, current_V, current_V + step_V, A*current_V+B, A*(current_V+step_V)+B);
						break;
					}
				}
				// special case, when HM jumped to PM
				if(myabs(PM_m - HM_m) < adapt_min_m_difference) { // we do have jump or II order transition
					if(myabs(step_V/adapt_step_reduction) >=  adapt_min_step) {
						step_V /= adapt_step_reduction;
						printf("#PM_m = HM_m -> smaller step! n=%f V=%f eF=%f step=%f: ", par.n(), par.V(), par.epsilonF(), step_V);
						printf("PM: E=%e, m=%e | HM: E=%e, m=%e, E_PM-E_HM=%e", res_PM.glob.Energy, res_PM.glob.magnetization, res_HM.glob.Energy, res_HM.glob.magnetization, res_PM.glob.Energy-res_HM.glob.magnetization);
						root_PM->print(); printf(" <PM|HM>  "); root_HM->print(); printf("\n"); fflush(stdout);
						continue;
					} else {
						printf("#");
						root_PM->print(); printf(" <PM|HM>  "); root_HM->print(); printf("\n"); fflush(stdout);
						printf("#PM_m = HM_m -> I got II order transition or jump between HM and PM\n");
						printf("#n=%f, V in (%f, %f); eF in (%f, %f)\n\n", current_n, current_V, current_V + step_V, A*current_V+B, A*(current_V+step_V)+B);
						break;
					}
				}

				// I'm checking if the step is too small
				//if(myabs(step_V) <  adapt_min_step) {
				//	printf("# BREAK Step is too small.\n\n");
				//	break
				//}

				// if you are here, it works, you can print stuff and as a starting point use the new variables
				current_V += step_V;
				printf("%f  %f  %f  %f", par.n(), par.V(), par.epsilonF(), step_V);
				root_PM->print();
				printf(" -2  ");
				root_HM->print(); printf("\n"); fflush(stdout);
				initV_PM.set(res_PM.get_V());
				initV_HM.set(res_HM.get_V());

				// decide if you should change the step
				int make_the_step_bigger = 0;
				int make_the_step_smaller = 0;
				double x, y;
				for(int k=0; k<6; k++) {
					//fork
					if(k==0) { x = res_PM.get_V()->mu(); y = old_PM_lv2.mu(); }
					else if(k==1) { x = res_HM.get_V()->mu(); y = old_HM_lv2.mu(); }
					else if(k==2) { x = PM_m; y = OLD_PM_m_lv2; }
					else if(k==3) { x = HM_m; y = OLD_HM_m_lv2; }
					else if(k==4) { x = PM_E; y = OLD_PM_E_lv2; }
					else if(k==5) { x = HM_E; y = OLD_HM_E_lv2; }

					if( myabs(x-y) < adapt_min_difference) { make_the_step_bigger++; }
					else if( myabs(x-y) > adapt_max_difference) { make_the_step_smaller++; }
				}

				if(make_the_step_smaller > 0) {
					step_V /= adapt_step_reduction;
				}
				else if(make_the_step_smaller == 0 && make_the_step_bigger > 0) {
					step_V *= adapt_step_reduction;
				}

				// copy new variables to OLD
				old_PM_lv2.set(res_PM.get_V());
				OLD_PM_E_lv2 = PM_E;
				OLD_PM_m_lv2 = PM_m;
				old_HM_lv2.set(res_HM.get_V());
				OLD_HM_E_lv2 = HM_E;
				OLD_HM_m_lv2 = HM_m;
		}
	}

	par.set_n(stop_n); par.set_V(stop_V); par.set_epsilonF(A*stop_V+B);
	printf("\n#END -- automatic_adaptive_run_nV\n");
	fflush(stdout);

	return;

}



} // end of namespace case2D2t_SGA_tPrim0
